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ABSTRACT 

This  report  summarises  recent  research  conducted  by  the  DSTO  in  plume 
model  development  for  urban  environments,  with  an  emphasis  on  establishing 
clear  physical  grounds  for  the  models,  yet  maintaining  enough  simplicity  to  be 
treated  numerically  in  an  operationally  viable  way.  The  aim  is  not  to  replace 
existing  operational  models  with  a  new  generation  of  more  accurate  models, 
but  to  provide  a  more  physics-based  framework  for  flow  and  dispersion  in  an 
urban  environment  that  can  reconcile  the  empirically  based  approach  of  current 
operational  models,  and  the  more  sophisticated  computational  fluid  dynamics 
techniques  now  gaining  popularity  for  atmospheric  dispersion  applications.  A 
key  feature  of  the  model  framework  developed  in  this  report  is  the  definition  of 
a  single  parameter  that  describes  canopy  morphology,  and  links  this  to  canopy 
flow  variables.  A  simple  canopy  dispersion  model  is  then  developed,  based  on 
flow  parameters  generated  by  the  canopy  model.  In  relevant  areas  the  well- 
known  Urban  Dispersion  Model  by  the  UK  Defence  Science  and  Technology 
Laboratory  is  used  as  a  benchmark  for  comparison.  Supporting  evidence  for 
the  models  developed  here  is  supplied  through  comparison  with  experimental 
data  from  a  water  channel  simulation. 
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Turbulent  Dispersion  Modelling  in  a  Complex  Urban 
Environment  —  Data  Analysis  and  Model  Development 

Executive  Summary 

The  purpose  of  this  report  is  to  summarise  the  recent  research  conducted  by  the  De¬ 
fence  Science  and  Technology  Organisation  (DSTO)  in  the  development  of  a  modelling 
framework  for  predicting  the  dispersion  of  chemical,  biological  or  radiological  (CBR)  con¬ 
taminants  in  urban  environments.  In  relevant  areas  the  well-known  empirically  based 
Urban  Dispersion  Model  (UDM)  developed  by  the  UK  Defence  Science  and  Technology 
Laboratory  (Dstl)  is  used  as  a  benchmark  for  comparison.  Comparisons  are  also  made 
with  some  more  practically  oriented  computational  fluid  dynamics  (CFD)  models  that  are 
receiving  considerable  interest  in  the  literature. 

The  aim  is  not  to  replace  existing  operational  models  with  a  new  generation  of  more 
accurate  models — the  broad  scope  of  issues  and  problems  to  be  addressed  in  completing 
such  a  task  is  not  able  to  be  attempted  here,  and  the  enhancement  of  operational  capability 
from  such  an  undertaking  is  far  from  clear.  Rather,  we  provide  a  more  physics-based 
framework  for  flow  and  dispersion  models  in  an  urban  environment  that  can  to  some 
extent  reconcile  the  very  empirically  based  approach  of  models  such  as  the  UDM,  and  the 
more  sophisticated  CFD  models.  In  developing  such  a  framework,  we  are  able  to  better 
interpret  experimental  data  for  validation  of  various  models,  and  provide  a  starting  point 
for  the  development  of  more  complex  concentration  fluctuation  models.  Such  models 
are  a  challenge  to  formulate  in  an  urban  environment,  but  ultimately  are  required  for 
the  development  of  synthetic  CBR  environments  to  predict  realistic  fluctuating  challenge 
levels,  toxicity  response,  optimal  detector  networks,  data  fusion  algorithms  for  enhanced 
CBR  situational  awareness,  and  the  rigorous  quantification  of  uncertainty  in  dispersion 
predictions. 

To  achieve  clarity  in  the  range  of  approaches  available  for  practical  models  of  con¬ 
taminant  dispersion  in  urban  areas,  an  overview  of  various  methods  recently  developed 
for  turbulent  transport  and  mixing  over  an  inhomogeneous  canopy  is  presented  in  Sec¬ 
tion  2.  The  analytical  complexity  of  those  methods  is  simplified  to  a  degree  that  allows 
straightforward  practical  implementation  and  application. 

Using  these  results  as  a  foundation,  in  Section  3  a  new  theoretical  framework  is  pro¬ 
posed  based  on  matching  flow  parameters  inside  and  outside  the  canopy.  The  proposed 
theoretical  framework  can  help  to  validate  and  justify  some  heuristic  assumptions  of  the 
UDM  and  may  be  used  as  a  valuable  performance  check  in  other  dispersion  models.  A 
key  feature  of  this  model  framework  is  the  definition  of  a  single  parameter  that  describes 
canopy  morphology,  and  links  this  to  canopy  flow  variables.  The  flow  model  is  used  as 
input  to  a  dispersion  model  for  the  canopy. 

In  Section  4  the  theoretical  findings  are  validated  against  experimental  data  collected 
from  water  channel  simulations  of  dispersion  through  a  variety  of  different  obstacle  arrays. 
These  arrays  provide  a  systematic  approach  to  investigating  the  effects  of  differing  obstacle 
morphology  on  flow  and  dispersion.  The  degree  to  which  the  model  can  reproduce  various 
plume  effects  due  to  the  canopy  morphology  is  discussed. 

The  framework  developed  allows  a  link  between  previously  published  DSTO  work  on 
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fluctuating  plume  models  in  urban  settings  [30,  33],  current  DSTO-Conrmonwealth  Scien¬ 
tific  and  Industrial  Research  Organisation  (CSIRO)  collaborative  work  on  internal  concen¬ 
tration  fluctuations  [37],  and  more  complex  distributed  drag  CFD  models  in  urban  areas 
[18,  19,  20].  This  will  enable  a  high  fidelity  model  of  urban  dispersion  to  be  implemented 
for  emergency  response  applications,  with  the  extra  capability  for  predicting  concentration 
fluctuations.  A  pilot  system  involving  coupling  between  a  nested  hierarchy  of  mesoscale 
meteorological  models,  CFD  models  and  Lagrangian  particle  models  has  already  been  es¬ 
tablished  by  Environment  Canada  in  collaboration  with  Defence  Research  &  Development 
Canada  (DRDC)  and  various  Canadian  universities  [2,  3].  The  methodology  developed 
here  serves  as  an  interim  step,  contributing  to  extra  capability  for  concentration  fluctua¬ 
tion  prediction  and  characterisation  of  the  uncertainty  in  modelling  predictions.  There  is 
also  potential  for  a  similar  system  to  be  implemented  within  Australia. 

The  framework  developed  has  also  provided  a  link  through  the  UDM  to  another  im¬ 
portant  operational  modelling  capability  that  DSTO  now  possesses — the  CBR  Virtual 
Battle-space  developed  by  Dstl.  There  is  further  scope  for  development  of  concentration 
realisation  models  within  this  capability  framework.  Finally,  the  generation  of  prototype 
simple  synthetic  environments  for  plume  realisations  is  also  required  to  further  develop 
data  fusion  algorithms  for  CBR  source  characterisation  and  network  detection.  The  frame¬ 
work  developed  here  is  the  first  step  in  the  development  of  such  a  synthetic  environment 
here  at  DSTO  for  these  purposes,  and  can  ultimately  be  used  to  compare  the  approaches 
studied  at  DSTO  with  those  of  our  international  collaborators. 
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1  Effect  of  a  canopy  on  the  wind  profile  and  the  resultant  model  of  the  mixing 

length  profile  lm.  The  wind  profile  above  the  canopy  is  approximately  the 
wind  profile  above  a  smooth  surface,  lifted  by  d,  as  in  (7) .  7 

2  Average  horizontal  velocity  profile.  The  pink/solid  line  is  for  a  sparse  canopy 

(e  =  0.5,  typical  urban)  and  the  blue/dashed  line  for  a  dense  canopy  (e  =  2, 
forest) .  14 


3  Left  -  Pink:  Relative  displacement  height  d/ H  as  a  function  of  canopy  drag 
parameter  e  (solid  pink  line).  That  d  <  0  for  e  — >  0  is  an  artifact  of  the 
functional  form  imposing  a  behaviour  change  in  flow  at  height  H,  even  for 
an  infinitely  sparse  canopy.  This  is  remedied  by  letting  H  — >  0  as  e  — - »  0. 

The  dashed  pink  line  is  a  sketch  of  d  corresponding  to  the  UDM  (11).  Left  - 
Blue:  Relative  roughness  height  zq  as  a  function  of  canopy  drag  parameter 
e.  Similarly  to  d,  the  fact  that  zo  ^  0  for  e  — ►  0  is  an  artifact  of  the  functional 
form  imposing  a  behaviour  change  in  flow  at  height  H,  even  for  an  infinitely 
sparse  canopy.  This  is  remedied  by  letting  H  — >  0  as  e  — >  0.  The  dashed 
blue  line  is  a  sketch  of  zq  corresponding  to  UDM  (12).  Right:  Plot  of  wind 
shear  stress  at  ground  level  computed  with  the  proposed  model.  We  see  a 
clear  monotonic  decrease  of  surface  shear  stress  with  the  canopy  density  e,  as 
intuitively  reasoned  on  physical  grounds .  15 

4  Plot  of  mixing  length  lm  for  our  flow  model.  The  pink  (long  dashed)  and  blue 
lines  are  for  e  =  1  and  4  respectively.  The  red  (short  dashed)  line  is  to  give  a 


comparison  of  the  gradient  of  Figure  21  in  [14]  for  z  >  H .  16 

5  Illustration  of  the  concept  of  a  virtual  source  induced  by  the  wake  regions  of 

the  obstacle  array.  The  virtual  source  is  lifted  to  height  CqH  (52) . 18 


6  Left:  Coanda  Urban  Array  001  consisting  of  regularly  spaced  cubic  objects 

with  packing  density  0.25.  Right:  Location  of  positions  within  a  cell  unit 
where  velocity  measurements  were  taken.  The  dark  shaded  area  indicates  the 
position  of  the  object  within  the  unit  cell . 20 

7  Vertical  cross-sections  of  the  normalised  mean  concentration  as  a  function  of 


downwind  position.  Solid  line  is  the  stretched  exponential  fit  (24).  Ca  depicts 
the  time  averaged  concentration  at  that  height  and  downstream  position, 
Ca,max  is  the  maximum  concentration  along  this  vertical  axis.  craz  is  the 
standard  deviation  in  the  vertical  position.  Description  of  the  obstacle  arrays 
can  be  found  in  [29] .  22 

8  Plot  of  log(log(Cmax/Ca))  against  normalised  vertical  coordinate  log  z/aaz- 

The  slopes  of  curves  are  determined  by  the  power  exponent  a  in  (24) . 23 

9  Variation  in  the  parameter  a  of  (24)  as  a  function  of  downstream  distance 

from  the  source  for  the  Coanda  data  with  no  obstacles .  24 
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10  Measurements  of  height  z  vs.  mean  stream- wise  velocity  U  for  the  Coanda 
Urban  Array  001.  The  data  points  come  from  measurements  at  different 
representative  points  around  a  single  array  object.  The  configuration  of  the 
array  and  the  location  of  measurement  points  can  be  seen  in  Figure  6  (right 
panel).  The  dashed  blue  line  indicates  the  height  of  the  objects  in  the  array. 

The  solid  line  is  a  fit  of  all  of  the  data  to  the  functional  form  given  in  (41) 

and  (42) . 24 

11  Plume  concentration  profiles  from  Coanda  data:  a  ground  level  release  over 

flat  terrain  (no  canopy) .  The  plots  on  the  left  show  the  vertical  concentration 
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Glossary 

[ABL]  Atmospheric  Boundary  Layer 
[CBR]  Chemical,  Biological  and  Radiological 
[CFD]  Computational  Fluid  Dynamics 

[CSIRO]  Commonwealth  Scientific  and  Industrial  Research  Organisation 

[DNS]  Direct  Numerical  Simulation 

[DRDC]  Defence  Research  and  Development  Canada 

[Dstl]  Defence  Science  and  Technology  Laboratory,  UK 

[NS]  Navier  Stokes 

[UDM]  Urban  Dispersion  Model 
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1  Introduction 

There  exists  a  large  variety  of  approaches  in  attempting  to  solve  the  difficult  problem  of 
atmospheric  dispersion  modelling  in  urban  environments.  Predicting  the  evolution  of  a 
plume  is  a  complex  problem,  to  be  described  by  advanced  methods  of  fluid  dynamics,  the 
theory  of  turbulence  and  diffusion,  and  statistics.  Its  comprehensive  modelling  is  compu¬ 
tationally  very  intensive  and  time  consuming.  However  in  an  operational  environment,  a 
prediction  is  needed  quickly. 

The  UK  Defence  Science  and  Technology  Laboratory  (Dstl)  Urban  Dispersion  Model 
(UDM)  [1]  provides  an  example  of  an  effective  operational  model  that  can  provide  quick 
answers  to  threat  scenarios  of  the  release  of  chemical,  biological  or  radiological  (CBR) 
material  of  gas  or  aerosol  form  into  an  urban  environment.  Understandably,  such  quick 
running  models  are  based  on  many  heuristic  and  practical  assumptions,  and  rely  heavily 
on  empirical  parameters  derived  under  certain  experimental  conditions  that  are  not  always 
clearly  extendable  to  a  wider  variety  of  atmospheric  conditions. 

On  the  other  hand,  there  has  been  growing  interest  in  the  use  of  computational  fluid 
dynamics  models  (CFD)  for  modelling  contaminant  dispersion  in  urban  areas.  Although 
more  physically  based,  these  models  tend  to  be  computationally  intensive,  and  require 
considerably  greater  computational  resources  (in  both  computing  power  and  time),  so 
that  they  are  not  generally  practical  for  operational  situations.  Despite  the  greater  em¬ 
phasis  on  the  underlying  physics  in  these  models,  the  inherent  complexity  of  flow  in  real 
urban  environments,  and  the  many  approximations  and  empiricism  still  required  in  the 
formulation  of  these  models  does  not  necessarily  imply  that  they  provide  significantly  more 
accurate  predictions. 

This  report  summarises  recent  work  at  DSTO  in  developing  a  modelling  framework 
to  describe  the  dispersion  of  CBR  contaminants  in  an  atmospheric  boundary  layer  (ABL) 
over  a  heterogeneous  surface.  The  aim  is  not  to  replace  existing  operational  models  with 
a  new  generation  of  more  accurate  models — the  broad  scope  of  issues  and  problems  to  be 
addressed  in  completing  such  a  task  is  not  able  to  be  attempted  here,  and  the  enhancement 
of  operational  capability  from  such  an  undertaking  is  far  from  clear.  Rather,  we  provide 
a  more  physics-based  framework  for  flow  and  dispersion  models  in  an  urban  environment 
that  can  reconcile  to  some  extent  the  very  empirically  based  approach  of  models  such 
as  the  UDM,  and  the  more  sophisticated  CFD  models  (being  developed  for  example  for 
Canadian  civilian  emergency  response  capability  [2,  3]).  In  developing  such  a  framework, 
we  are  able  to  better  interpret  experimental  data  for  validation  of  various  models,  and 
provide  a  starting  point  for  the  development  of  more  complex  concentration  fluctuation 
models.  Such  models  are  a  challenge  to  formulate  in  an  urban  environment,  but  ultimately 
are  required  for  the  development  of  synthetic  CBR  environments  to  predict  realistic  fluc¬ 
tuating  challenge  levels.  Concentration  fluctuations  are  important  in  analysing  toxicity 
response,  optimal  detector  networks,  data  fusion  algorithms  for  enhanced  CBR  situational 
awareness,  and  the  rigorous  quantification  of  uncertainty  in  dispersion  predictions. 

In  order  to  give  a  background  understanding  of  previous  work  in  this  area,  in  Section  2 
a  brief  overview  of  the  complex  modelling  methods  recently  developed  for  turbulent  trans¬ 
port  and  mixing  is  presented.  In  Section  3  a  new  theoretical  framework  is  proposed  based 
on  matching  flow  parameters  inside  and  outside  the  canopy.  The  proposed  theoretical 
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framework  can  help  to  validate  and  justify  some  heuristic  assumptions  of  the  UDM,  and 
may  be  used  as  a  valuable  performance  check  in  other  dispersion  models.  In  Section  4  data 
from  our  water  channel  experimental  data  are  used  in  validating  our  theoretical  findings. 


2  An  Overview  of  Turbulent  Dispersion  in  the 

Urban  Boundary  Layer 


2.1  Velocity  Field  in  the  Urban  Boundary  Layer 


Since  flow  in  the  ABL  is  a  key  driver  for  CBR  pollutant  advection  and  dispersion,  the 
development  of  a  high  fidelity  model  of  this  flow  is  a  crucial  step  in  the  modelling  of  the 
whole  turbulent  dispersion  process. 

The  simplest  model  of  flow  in  the  ABL  is  the  celebrated  log-law  profile  [4].  With 
coordinate  x  the  downstream  distance  from  the  source,  y  is  the  crosswind  coordinate,  and 
z  >  0  the  height  above  the  ground,  the  vertical  profile  of  mean  flow  aligned  with  the  x-axis 
takes  the  form 

c,(2)  =  ^loe(^)'  (1) 

where  tt*  is  the  friction  velocity,  k  =  0.41  is  von  Karman’s  constant,  zo  =  v/u*  is  the 
viscous  length1,  and  v  is  the  viscosity.  Equation  (1)  is  valid  in  the  surface  layer  portion 
of  the  ABL,  for  neutral  stability  conditions. 

It  has  been  known  for  a  long  time  (dating  back  to  Prandtl  [4])  that  the  ABL  velocity 
profile  can  be  fairly  approximated  by  a  power- law  function: 

{  z\m 

U(z)  =  au*y—J  ,  (2) 


where  a  and  m  are  constants.  This  profile,  being  algebraically  simpler  than  the  log-law 
profile,  has  been  used  merely  as  a  convenient  engineering  approximation,  but  recently  it 
has  attracted  much  attention  after  being  shown  that  it  can  be  justified  based  on  a  self¬ 
similarity  property  of  the  ABL  flow.  For  the  boundary  layer  over  a  flat  smooth  surface,  it 
was  rigorously  shown  by  Barenblatt  et  al  [5,  6]  that 


log  Re  5  3 

=  U7U  +  2’  m  =  UZr? 


(3) 


where  Re  is  the  Reynolds  number  of  the  flow.  For  our  purposes  we  assume  that  both 
profiles  (1)  and  (2)  hold  equally  for  turbulent  flow  in  the  ABL  with  parameters  a  and 
m  to  be  determined  from  data  fitting.  This  is  justified  due  to  the  stochastic  nature 
of  atmospheric  flows.  The  flow  variables  in  the  atmosphere  do  not  remain  statistically 
stationary  for  long  enough  for  truly  converged  mean  quantities,  so  the  resulting  uncertainty 
in  measured  profiles  of  mean  velocity  usually  cannot  provide  a  distinction  between  (1)  and 
(2). 

1for  flow  over  smooth  surfaces,  or  roughness  length  for  flow  over  rough  surfaces 
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Observed  values  of  m  in  the  atmosphere  range  from  nearly  0  in  very  unstable  condi¬ 
tions,  representing  perfect  mixing  and  a  uniform  velocity  profile,  to  nearly  1  in  very  stable 
conditions,  approaching  the  Couette  linear  profile  of  laminar  motion  over  a  plane  surface. 
In  neutral  conditions  m  ~  1/7  [4]  (a  thorough  discussion  of  wind  profiles  in  different 
weather  conditions  can  be  found  in  [35]).  The  value  of  m  also  depends  on  surface  rough¬ 
ness:  roughness  promotes  mixing  near  the  surface,  which  reduces  the  velocity  gradient  at 
small  z  and  thus  leads  to  larger  m.  These  effects  can  be  computed  in  a  semi-empirical  way 
and  the  UDM  [1]  also  provides  a  framework  for  this.  While  it  is  acknowledged  that  the 
two  profile  forms  are  not  mathematically  equivalent,  we  use  both  forms  in  various  parts  of 
this  report  to  derive  tractable  analytical  models  of  flow  and  dispersion  within  urban  areas 
that  would  not  be  possible  by  restricting  ourselves  to  just  one  choice. 

The  turbulent  (or  eddy)  viscosity  in  the  ABL  can  be  determined  as 

^) =*»(£)”■  « 

where  Kq  and  r  are  constants.  The  value  of  n  can  be  specified  based  on  various  analytical 
models  of  turbulent  boundary  layers  or  directly  estimated  from  the  observational  data.  In 
particular,  for  so-called  conjugate  power-law  profiles  (constant  momentum  flux)  [4] 

n  =  1  —  m,  (5) 

and  for  Monin-Obukhov  similarity  profiles 

n  =  1.  (6) 

In  general,  we  may  still  consider  n  as  a  free  independent  parameter  of  the  flow  model. 

Over  a  rough  surface  (i.e.  ABL  over  canopy)  the  profiles  (1)  and  (2)  are  modified 
according  to 


where  now  zq  is  the  roughness  height  and  d  is  the  so-called  displacement  height.  Both  d 
and  zq  should  be  considered  as  new  fitting  parameters  of  the  ABL  flow  over  the  canopy, 
and  their  simple  physical  interpretation  is  not  always  obvious.  For  instance  (7)  still  pro¬ 
duces  meaningful  results  even  for  d  <  0,  when  its  naive  interpretation  as  the  “ABL  flow 
displacement  by  canopy”  is  not  feasible.  It  is  worth  pointing  out  that  determining  zo  and 
d  simultaneously  from  experimental  data  is  problematic.  The  difficulty  arises  from  the  log 
relationship,  which  causes  an  exponential  dependency  of  the  parameters  on  uncertainties 
in  the  measured  data.  This  problem  has  been  described  in  more  detail  by  Gailis  [7],  and 
demonstrated  rigorously  in  a  Bayesian  inference  framework  for  certain  data-sets  [8]. 

The  ongoing  challenge  for  any  plausible  model  of  ABL  flow  over  the  heterogeneous 
canopy  is  to  establish  reliable  functional  relationships  between  zo,  d  and  the  parameters 
of  the  canopy  (i.e.  building  morphology).  Once  such  relationships  are  found,  they  can 
be  used  in  conjunction  with  (7)  or  (8)  to  model  the  mean  velocity  profile  affected  by  a 
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particular  canopy.  Such  a  modified  velocity  profile  can  later  be  used  as  an  input  for  the 
pollutant  transport  equation  (21). 

Various  simple  approaches  based  on  empirical  relations  can  be  found  in  the  literature 
to  relate  canopy  parameters  to  velocity  profile  parameters.  These  approaches  tend  to 
use  formulas  with  some  physical  basis,  but  often  coupled  with  heuristic  assumptions  and 
heavily  reliant  on  tuning  the  relationships  to  specific  data.  Two  older  examples  include 
those  of  Lettau  [9]  and  Counihan  [10].  The  approach  used  in  the  UDM  [1]  was  developed 
by  Macdonald  et  al  [11],  and  attempts  to  base  the  relationships  on  firmer  physical  grounds. 
The  first  key  canopy  parameter  is  the  ratio  of  the  collective  frontal  area  of  the  obstacles 
to  the  collective  horizontal  area  of  the  canopy  (called  the  lot  area): 


Xf  ~ 


_  f4/ 

AT 


(9) 


i.e.  the  density  of  the  windward  faces  of  the  canopy  (which  is  equivalent  to  leaf  area 
density  for  plant  canopies).  The  second  key  canopy  parameter  is  the  ratio  of  the  summed 
horizontal  area  of  objects  in  the  canopy  (plan  area)  to  the  lot  area 

A  =  (10) 


or  the  dimensionless  plan  area.  Then,  given  the  averaged  canopy  height  H,  and  the  drag 
coefficient  Cd  of  individual  canopy  elements,  some  simple  considerations  of  shear  stress 
give  the  following  relationships  [11]: 


d_ 

H 

zo 

H 


1  +  A~xp(Xp  -  1), 


(11) 

(12) 


Here  k  =  0.41  is  von  Karman’s  constant,  and  we  have  two  constants  A  and  f3  that  are 
tunable  according  to  particular  obstacle  arrays.  For  a  regular  cubic  array,  A  =  3.59  and 
(3  =  0.55,  whereas  for  a  staggered  cubic  array,  A  =  4.43  and  /?  =  1.0.  We  observe  that 
d  =  0  and  zq  =  0  for  Xp  =  0  (no  canopy) ,  and  d  =  H  and  zo  =  0  for  Xp  =  1  (full  packing) . 

The  formulae  (8),  (11)  and  (12)  are  the  core  components  of  the  surface  layer  flow 
model  in  the  UDM  and  plots  for  d  and  zq  corresponding  to  these  formulae  are  presented 
in  Figure  25  of  the  technical  documentation  [1].  In  Section  3.1  we  will  provide  a  theoretical 
framework  in  order  to  derive  some  analytical  relationships  between  zq,  d  and  the  parame¬ 
ters  of  the  canopy.  This  will  help  to  validate  and  justify  some  of  the  heuristic  assumptions 
and  tuning  used  to  generate  the  relationships  (11)  and  (12)  as  employed  in  the  UDM,  and 
may  be  used  as  a  valuable  performance  check  in  other  dispersion  models. 

The  approach  outlined  above  is  simple  and  uses  some  physics-based  arguments,  but 
is  primarily  given  its  predictive  power  through  empirical  relationships.  More  advanced 
approaches  rely  on  utilising  the  full  fluid  dynamics  equations  of  the  canopy  flows — usually 
referred  to  as  the  Navier  Stokes  (NS)  equations — with  various  simplifying  approximations 
and  empirical  inputs.  The  entire  range  of  techniques  in  computational  fluid  dynamics 
(CFD)  is  very  broad  and  beyond  the  scope  of  this  report,  but  here  we  focus  on  one 
approach  that  has  already  been  demonstrated  as  useful  for  practical  applications  involving 
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relatively  stringent  time  and  computational  requirements.  This  more  advanced  approach 
to  surface  layer  flow  in  the  canopy  is  to  compute  the  mean  flow  by  solving  the  time  and 
space  averaged  momentum  equations.  Thus  influence  of  the  canopy  is  taken  into  account 
in  a  bulk  fashion  (as  opposed  to  detailed  point  by  point  modelling)  through  its  viscous 
drag  and  form  (pressure  difference)  drag  on  the  flow,  which  consequently  has  led  to  the 
name  “distributed  drag”  for  this  class  of  models. 

Relatively  straightforward  implementations  of  the  distributed  drag  approach,  includ¬ 
ing  a  systematic  description  of  various  flow  regions/features  in  a  canopy,  and  physical 
approximations  that  may  be  made  within  these  regions  have  been  given  by  Belcher  et 
al  [13]  and  Coceal  et  al  [14,  15].  A  more  complex  approach  is  to  implement  a  modi¬ 
fied  spatially  averaged  k-e  CFD  model,  which  includes  budget  equations  for  the  turbulent 
kinetic  energy  k  and  turbulent  dissipation  e,  as  demonstrated  by  Katul  et  al  [16].  The 
above  mentioned  studies  have  all  made  some  implicit  assumptions  in  the  derivation  of 
the  spatially  and  Reynolds  averaged2  NS  equations,  most  importantly  a  neglect  of  the 
dispersive  stresses3  arising  from  the  spatial  averaging  operation.  Lien  et  al  [18,  19]  put 
the  mathematical  framework  in  dealing  with  time  and  spatial  averaging  of  NS  and  k-e 
equations  on  a  rigorous  footing,  and  demonstrated  the  non-commutation  of  the  spatial 
and  time  averaging  operators,  which  led  to  different  structural  equations  ultimately  de¬ 
scribing  the  same  physics,  depending  on  whether  spatial  or  time  averaging  is  applied  first. 
By  comparing  results  with  high  resolution  Reynolds-averaged  k-e  NS  models,  they  also 
demonstrated  that  the  dispersive  stresses  were  at  times  of  great  importance  in  canopies 
with  urban  morphology4  [20] . 

For  the  purposes  of  this  report,  it  is  sufficient  to  take  the  following  form  for  the  mean 
momentum  equations,  without  elaborating  on  the  mathematical  details  for  the  issues  raised 
in  the  previous  paragraph: 


dvi  dvi  _  dp  dTij  |v|uj 

-'I.  “  ^  W  .'1  .'1  -I  ./  i  ■  ./  i  r 

at  oxj  oxi  oxj  Lc 

where  Tij  is  a  spatially- averaged  Reynolds  stress,  and  ft  is  a  smoothly- varying  canopy 
element  drag  (which  arises  from  spatially  averaging  the  localised  drag  due  to  individual 
roughness  elements).  The  parameter 


2  H  (1-Ap) 
CD(z)  A  f 


(14) 


does  not  depend  on  velocity  and  is  called  the  canopy  drag  length  scale;  as  before,  H  is 
the  average  height  of  canopy  objects,  and  A j  and  Xp  are  morphological  parameters  of  the 
canopy.  The  physical  meaning  of  Lc  is  a  length-scale  for  an  incident  wind  profile  to  adjust 
to  the  canopy  (for  typical  values  of  Lc,  see  [13]). 

2  Reynolds  averaging  is  ensemble-  averaging  and  reduces  to  time  averaging  for  statistically  stationary 
flows. 

3  Specifically,  the  dispersive  stresses  arise  from  the  spatial  correlation  in  the  time-mean  velocity  field  as 
it  varies  with  position  in  an  averaging  volume. 

4Other  studies  of  urban  canopies  have  asserted  that  dispersive  stresses  are  negligible  by  quoting  data 
in  the  literature  taken  from  plant  canopies.  Although  this  may  be  correct  for  plant  canopies,  it  will  be 
demonstrated  below  (and  similarly  has  been  stated  in  the  urban  canopy  literature)  that  morphological 
parameters  are  considerably  different  between  urban  and  plant  canopies,  so  that  conclusions  drawn  from 
plant  canopy  data  are  not  necessarily  transferable. 
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The  formula  (14)  often  permits  further  simplification  by  specifying  height-averaged 
Cd-5  For  the  flow  over  a  canopy  consisting  of  a  cubical  array  Cd  ~  2  [13]  (Cd  ~  1.2  in 
the  UDM  [l]6 *),  so  (14)  can  be  reduced  to 

(15) 

where  we  introduced  the  non-dimensional  “canopy  density”  parameter  e  for  convenience. 
Since  e  is  the  only  parameter  determined  by  the  canopy  in  the  “distributed  drag”  frame¬ 
work,  it  is  the  parameter  that  should  be  used  for  any  canopy  parameterisation.  In  Sec¬ 
tion  3  we  will  analytically  derive  expressions  for  d  and  zq  via  e. 

It  is  evident  that  canopies  with  higher  e  are  denser  (for  example,  the  case  of  Xp  =  1 
corresponds  to  e  =  oo  i.e.  when  the  whole  underlying  surface  is  occupied  by  the  canopy). 
In  general  e  ~  1  for  urban  canopies  and  e  >  1  for  plant  canopies.  In  sparse  canopies, 
horizontal  momentum  flux  from  the  ABL  flow  is  not  influenced  much  by  the  canopy, 
and  it  directly  applies  to  the  underlying  surface;  otherwise  in  dense  canopies  all  drag  is 
fully  absorbed  by  the  canopy  elements.  In  some  respects  a  sparse  canopy  case  is  more 
challenging  since  it  requires  imposing  the  correct  boundary  condition  for  mean  velocity 
U(z)  on  the  underlying  surface.  This  depends  of  the  fraction  of  drag  propagating  down  to 
the  surface  and  is  not  known  beforehand  (in  dense  canopies  we  simply  have  dU(z)/dz  =  0 
at  z  =  0). 

To  illustrate  the  significance  of  the  canopy  parameters  described  in  this  framework 
with  a  simple  example,  we  consider  a  canopy  consisting  of  a  regular  array  of  identical 
cylinders  with  radius  ro  and  height  H,  separated  by  a  distance  of  2ro-8  In  this  example, 
the  expression  (15)  undergoes  further  simplification,  with  A /  =  ( 2H  /irr^Xp ,  and  Cd  ~  1, 
so  the  formula  can  be  reduced  to 

H  Xp 

e  ~ - - — 

2irro  1  —  Xp 

We  can  see  how  the  value  of  parameter  e  depends  on  the  contribution  of  two  factors: 
the  shape  of  individual  elements  of  the  canopy  and  their  packing  density.  None  of  these 
parameters  separately  determines  the  limiting  value  of  e  in  the  case  Xp  — >  1  or  H  — * 
0,  r o  — ►  0.  In  Sections  3.1  and  Appendix  A  we  argue  that  the  limit  should  be  taken  in 
conjunction  with  H  —>  0.  Of  course,  this  formula  can  be  used  for  approximate  estimation 
of  e  for  a  canopy  consisting  of  other  elements,  not  just  cylinders. 

In  order  to  derive  any  analytical  solution  of  (13),  a  closure  assumption  must  be  used. 
The  spatially-averaged  Reynolds  stress  is  usually  represented  with  eddy  viscosity  models 

5This  is  particularly  important  in  plant  canopies,  where  considerably  different  drag  is  experienced  at 
different  heights  (the  crowns  of  trees  versus  trunks),  but  is  less  significant  for  the  mainly  uniform  obstacle 
arrays  considered  in  experimental  studies  of  urban  dispersion,  such  as  considered  in  this  paper. 

6This  highlights  a  particular  advantage  of  the  distributed  drag  approach,  whereby  variation  in  canopy 
density  and  consequently  variation  in  canopy  drag  can  be  explicitly  incorporated  in  the  formalism,  whereas 
for  simple  empirical  models  such  as  the  UDM,  a  particular  value  had  to  be  selected  a  priori. 

'While  the  equations  for  e,  such  as  (15),  and  that  derived  in  section  3.1,  produce  broadly  correct 
canopy  drag  behaviour,  they  do  not  capture  all  subtleties  created  by  non-linear  effects  of  currents  between 
buildings.  This  is  a  case  of  a  practical  approximation  being  made  in  order  to  generate  results  that  can  be 
tractably  used  and  that  are  relatively  universal  compared  to  CFD  calculations  of  the  drag. 

8  This  array  is  like  the  regular  cubical  obstacle  array  discussed  in  the  Coanda  data-set  (see  Section  4), 
but  with  the  cubes  replaced  by  cylinders. 
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Figure  1:  Effect  of  a  canopy  on  the  wind  profile  and  the  resultant  model  of  the  mixing 
length  profile  lm.  The  wind  profile  above  the  canopy  is  approximately  the  wind  profile  above 
a  smooth  surface,  lifted  by  d,  as  in  (7). 


Tij  K(dvi/dxj),  and  variations  of  K-e  turbulence  models  can  be  employed  for  closure 
(for  details  see  [13,  14,  15,  16,  17,  18,  19,  20]).  Among  the  variety  of  available  closures,  the 
mixing  length  hypothesis  has  been  widely  adopted  on  the  grounds  of  its  simplicity  even 
though  “the  assumptions  behind  this  hypothesis  are  rarely  met”  [17]. 

In  the  mixing  length  approach,  the  whole  complexity  of  the  turbulent  flow  parameter- 
isation  is  reduced  to  a  heuristic  model  for  the  mixing  length  function  lm(z).  The  mixing 
layer  approach  has  some  justification  in  the  case  of  a  boundary  layer  over  a  flat  smooth 
plane  where  it  produces  a  very  close  match  to  experimental  data,  but  it  is  evident  that  by 
its  definition  the  mixing  length  profile  is  not  a  generic  property  and  it  may  be  very  canopy 
specific.  It  is  also  very  different  for  sparse  and  dense  canopies:  in  a  dense  canopy  the 
turbulent  flow  in  the  ABL  cannot  penetrate  deeply  into  the  canopy  and  affect  the  mixing 
process  near  the  underlying  surface,  while  in  a  sparse  canopy  mixing  near  the  surface  is 
similar  to  that  of  a  free  surface.  Below  are  some  examples  of  the  models  for  lm  employed 
in  [13,  14,  15,  16].  The  profile 


1  1  1  , 

-  =  —  +  lc  =  const 

Lm  &  iq 


(17) 


was  used  for  an  urban  area  where  the  K  used  in  the  literature  is  a  non-dimensional 
constant,  and 


lm  =  [Z-  (2/3 )H)K,  z/H  >  1,  (18) 

lm  =  const,  z/H  <  1  (19) 

for  plant  canopies  (see  Figure  1). 

An  interesting  result  was  recently  report  by  Coceal  et  al  [14],  where  the  function  of 
lm(z)  was  a  posteriori  evaluated  from  the  canopy  flow  calculated  by  Direct  Numerical 
Simulation  (DNS — where  all  scales  of  turbulence  are  explicitly  resolved  in  the  numerical 
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solution  of  the  NS  equation).  The  linear  sections  of  the  function  lm(z)  were  evident  in 
their  results  (see  Figure  21  of  [14]),  it  is  claimed  that  in  these  regions  the  boundary  layer 
profile  can  be  represented  in  the  “displacement”  form  (7).  Various  approximations  of  the 
lm  profile  presented  in  Figure  1  can  be  used  for  more  plausible  models  of  the  mean  flow 
in  a  canopy  array  in  the  mixing  length  parameterisation  approach.  In  dense  canopies  the 
variability  of  the  bottom  part  of  the  profile  disappears  since  here  lm  is  constant  for  e<l 
(see  (19)).  In  Section  3,  the  mixing  layer  profile  lm  based  on  our  model  of  the  mean  flow 
in  the  canopy  is  calculated. 

Invoking  the  mixing  length  hypothesis  with  lm  =  const  allows  an  analytical  solution 
for  the  velocity  profile  near  the  top  of  the  dense  canopy  (f>l)  that  is  known  as  Cionco’s 
canopy  profile  [21]:  U(z)  ~  exp  ( cz )  with  c  a  constant.  This  profile  is  used  in  the  UDM  as  a 
mean  flow  profile  in  the  canopy.  It  is  worth  noting  that  while  Cionco’s  profile  describes  well 
the  flow  near  the  top  of  the  canopy,  it  violates  the  boundary  condition  on  the  underlying 
surface  (1/(0)  /  0),  and  this  can  be  justified  only  for  a  dense  canopy. 

Another  approach  to  deriving  analytical  solutions  of  (13)  is  to  employ  various  hy¬ 
potheses  of  the  functional  relation  between  K  and  U  in  the  canopy.  Such  hypotheses  are 
usually  more  generic  than  the  simple  mixing  length  approach  since  they  are,  by  definition, 
universal  and  non-local  (contain  no  explicit  dependencies  of  space  coordinates  or  canopy 
shape).  For  instance,  the  closure  K  ~  U  leads  to  Cowan’s  profile:  U(z)  ~  y/smh(A~?) 
[38]  with  A  a  constant.  The  Cowan’s  velocity  profile  meets  the  boundary  condition  at  the 
bottom  of  the  canopy  and  we  will  use  Cowan’s  profile  as  an  universal  flow  model  in  the 
canopy.  In  Section  3.1  we  extend  Cowan’s  profile  for  the  mean  velocity  profile  in  the  ABL 
for  a  linear  K-U  relationship  that  is  valid  for  dense  and  sparse  canopies  (i.e.  for  large  and 
small  values  of  e),  determining  the  constants  of  the  model  in  terms  of  the  canopy  density. 

The  important  aspect  of  realistic  models  of  turbulent  diffusion  in  the  ABL  is  to  cater 
for  the  ABL  stability  effects.  These  effects  are  usually  incorporated  in  dispersion  models 
based  on  the  Monin-Obukhov  similarity  theory.  According  to  this  theory,  the  stability 
effects  can  be  accounted  for  by  alteration  of  mean  velocity  profiles  in  the  following  way: 

u(z)  =  U(z)  -  (t>m  ’  (20) 

where  L  is  the  so-called  Monin-Obukhov  length  scale  (i.e.  the  scale  of  turbulent  flow  with 
thermal  fluxes),  4>m(z)  is  a  known  dimensionless  function  [4,  17],  and  U{z)  is  either  of  the 
profiles  (1)  or  (2).  This  considerably  complicates  the  ABL  flow  problem,  and  will  not  be 
further  considered  in  the  current  report,  but  is  a  challenging  extension  for  a  more  general 
model  than  that  discussed  in  Section  3.1,  and  will  be  developed  in  future  work.  Progress 
for  plant  canopy  models  has  also  been  recently  made  by  Harman  and  Finnigan  [17]. 


2.2  Mean  Concentration  Field 

The  equation  for  mean  concentration  C  is  a  well-known  equation  of  the  pollutant  budget 
in  K-theory  [4,  22]: 


(21) 
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with  the  assumption  that 


K, 


yy(z)  —  2 


do|(a 

dx 


and  <Jy(x)  is  the  horizontal  plume  spread.  It  is  known  [4,  5,  6,  22,  23]  that  for  power-law 
profiles  and  a  source  on  the  ground  that  this  equation  provides  the  analytical  solution  for 
the  pollutant  concentration  in  the  surface  layer.  The  solution  can  be  written  in  the  form 


C(x,y,z)  =  Cy(x,y)Cz(x,z).  (22) 

The  horizontal  transverse  mean  concentration  profile  Cy  has  a  Gaussian  solution 

c'»(x’!,)  =  ^b;exp(-||)’  (23) 

where  ay  =  ay{x)  is  the  lateral  dispersion  parameter.  The  vertical  mean  concentration 
profile  Cz  has  a  stretched  exponential  solution 


Cz  =  C0(x)  exp  (-B(x)C)  ■ 

(24) 

(25) 

zq  x  azKo  a 

(26) 

and  a  =  2  +  m  —  n,  uq  =  au *,  T(/3)  is  the  Gamma  function,  and  Q  is  the  rate  at  which 
the  source  releases  the  pollutant  (rate  of  injection).  The  analytical  solution  (22),  (23)  and 
(24)  is  a  foundation  for  our  model.  In  Section  4.5  we  comprehensively  validate  it  against 
our  experimental  data. 

For  conjugate  (5)  or  Monin-Obukhov  (6)  profiles  we  have  simple  estimates 

ac  =  1  +  2m,  (27) 

cxmo  =  1  +  m.  (28) 

For  small  a  (neutral  conditions)  these  estimates  are  very  close,  but  in  general  the  following 
inequality  seems  to  hold:  aMO  <  a  <  olc.  This  estimate  will  be  validated  based  on  our 
experimental  data-set  (see  Section  4.4). 

One  of  the  important  properties  of  pollutant  dispersion  in  the  ABL  described  by  solu¬ 
tion  (22)  is  its  universal  scaling,  i.e.  the  simple  functional  relations  between  its  parameters: 

crz  ~  <Jy  ~  Ah  (29) 

where  y  is  the  plume  centroid  position  in  the  vertical  direction  (the  centre  of  gravity  of 
the  concentration  distribution).  Such  scaling  can  be  justified  based  on  properties  of  flow 
in  the  ABL  [22]  and  can  significantly  simplify  expressions  (23)  and  (24),  reducing  them  to 
a  self-similar  universal  form  in  dimensionless  variables. 

The  solution  (22),  (23)  and  (24)  can  be  generalised  for  an  elevated  source  of  height 
z  =  Hs  [4,  22,  23].  In  such  cases  Cy(x,y )  will  not  change,  but  Cz(x,z )  can  be  written  in 
the  form 

C,  =  Co(x,C)exp(-B(x)(C0“  +  r))/-,(2S(x)(CCo)a/2),  (30) 
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C0(x,0  =  —((tor/2aB(x), 


UqZq 


where 


/- /  \  z  Hs  x o  u0z$  n  (1  +  m) 

C(")  =  j  Co  =  — ,  B(x)  =  —,  x0  =  ^FF,  v=\- 


(31) 


(32) 


zq  zq  '  x  a2Ko  a 

a  =  2  +  m  —  n  as  above,  and  Iv{z )  is  the  type  I  modified  Bessel  function.  It  is  worth 
noting  that  parameter  £o  (source  elevation)  is  responsible  for  the  vertical  concentration 
profile  variability  close  to  the  source,  and  its  contribution  diminishes  at  x  >  Hs.  It  can 
be  shown  that  at  the  limit  Hs  — >  0  (i.e.  Co  — ►  0)  this  solution  tends  to  (24). 

For  the  solutions  given  by  (24)  or  (30),  the  pollutant  concentration  on  the  surface 
(C  =  0)  is  given  by 

n 

(33) 


C(Z)  =  -exp{-F/Q, 

with  £  =  x/xq  and  D,  F  constant.  We  can  see  that  function  C(£)  increases  from  zero, 
reaches  its  maximum,  and  then  decays  as  a  power-law.  The  exponent  7  of  power  law 
decay  is  determined  by  the  velocity  profile  (i.e.  exponent  m).  For  the  conjugate  profile 
the  exponent  7  is  given  by 

1  +  m  , 

7  =  - .  (34) 

[  1  +  2m  v  ’ 

It  should  be  noted  again  that  the  value  of  m  and  hence  7  is  also  influenced  by  the  type  of 

canopy  and  the  stability  conditions. 


In  the  UDM  framework,  a  Gaussian  form  is  used  for  the  concentration  profile,  i.e. 
simply  a  =  2  in  (24)  or  (30),  and  the  surface  reflection  is  included  for  elevated  sources  [1]. 
This  simply  leads  to  7  =  2  in  (33)  and  any  dependency  of  7  on  wind  shear  disappears. 
This  significantly  deviates  from  the  stretched  exponential  model  (24)  and  (30),  even  in  the 
neutral  case,  since  by  (34)  7  ~  1  for  m  <C  1.  This  may  result  in  possible  underestimation 
of  downstream  concentration  levels  predicted  by  UDM,  mostly  at  far  distances  from  the 
source.  This  deviation  will  only  increase  if  we  take  into  account  stability  conditions  and 
canopy  effects,  since  they  correspond  to  the  higher  values  of  m  causing  7  to  further  decrease 
and  to  approach  its  limiting  value  7  ~  1/2  following  from  (34).  Hence  the  effect  of  wind 
shear  may  be  very  important  for  CBR  concentration  prediction  and  should  be  taken  into 
account  for  high  fidelity  dispersion  modelling.9 


3  Development  of  a  New  Plume  Model 

3.1  Flow  in  the  Canopy 

In  Section  2.1  the  atmospheric  velocity  and  diffusivity  profiles  within  the  ABL  were  dis¬ 
cussed  for  the  case  with  very  little  or  no  canopy  to  impede  the  flow  based  on  (1)  and  (2). 

9This  effect  may  not  be  noticeable  in  the  urban  modules  of  the  Hazard  Prediction  and  Assessment 
Capability  (HPAC)  of  which  UDM  forms  a  component,  as  when  the  plume  reaches  about  1  km  or  10-12 
obstacle  rows  downwind  of  the  source,  UDM  passes  puffs  off  to  the  far  field  dispersion  model  SCIPUFF. 
Another  complication  in  this  comparison  is  that  UDM  is  a  puff  model,  and  thus  though  the  vertical  profile 
of  puffs  may  not  be  optimal  according  to  the  gradient-diffusion  theory  described  here,  an  aggregation  of 
hundreds  or  thousands  of  puffs  in  a  computational  domain  may  not  lead  to  the  simple  £~2  dependency  for 
C(£)  given  by  the  plume  framework  of  (33). 
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The  case  of  wind  flow  above  a  significant  canopy  was  also  discussed,  viz.  (7)  and  (8). 

Here  the  canopy  wind  flow  model  of  Cowan  [38]  is  extended  to  create  a  self-consistent 
model  for  the  full  profile,  both  within  and  above  the  canopy.  The  aim  in  the  development 
of  this  model  is  to  employ  as  few  heuristic  assumptions  as  possible,  and  for  the  model 
to  be  dependent  on  as  few  external  parameters  as  is  physically  reasonable.  Supporting 
experimental  evidence  is  given  in  Section  4.4. 

The  key  to  modelling  average  velocity  flow  within  a  canopy  is  to  study  the  shear  stress 
within  the  canopy,  caused  by  the  drag  the  canopy  exerts  on  the  flow  [38,  24,  25].  In  this 
way,  we  can  derive  a  common  framework  to  link  equations  employed  in  empirically  based 
and  heuristic  models  such  as  the  UDM,  and  more  physics  based  distributed  drag  CFD 
models  described  in  Section  2.1.  The  shear  stress  is  given  by 

t  =  PK(z)dU^\  (35) 

where  p  is  the  fluid  density,  and  I\(z )  =  I\zz  and  U  are  the  average  velocity  and  diffu- 
sivity  (as  in  Section  2.1).  It  is  important  to  note  that  in  the  K-theory  of  dispersion,  the 
assumption  of  equality  between  eddy  viscosity  and  eddy  diffusivity,  namely  K(z )  =  Kzz , 
is  automatically  made,  but  this  is  by  no  means  grounded  in  rigorous  physical  theory.  We 
will  follow  the  same  approach  in  this  section  to  enable  the  formulation  of  tractable  ana¬ 
lytic  models  with  minimal  fitting  parameters,  however  in  general  it  may  be  assumed  that 
a  constant  of  proportionality  other  than  unity  exists. 

The  rate  of  change  of  shear  stress  through  the  canopy  can  be  related  to  the  drag  by 

^  =  PCdAU\z ),  (36) 

where 

A=^~  =  X±- 
AtH  H 

is  the  area  density  (the  amount  of  area  perpendicular  to  the  flow  per  unit  volume).  It  is 
assumed  that  the  area  density  is  constant  in  z,  a  reasonable  approximation  for  an  urban 
canopy.  Assuming  self-similarity  between  U(z)  and  K(z)  within  the  canopy10  [25]  gives 

K(z)  =  aHU(z)  (38) 


(37) 


for  canopy  height  H ,  with  a  a  dimensionless  parameter  to  be  determined  by  the  model 
(see  below).  This  leads  to  the  following  differential  equation  describing  how  drag  affects 
velocity  within  the  canopy: 


2e 


+  2  U(z) 


d  2U(z) 
d  z'2 


(39) 


Here  we  see  the  re-emergence  of  the  parameter  e  defined  in  (15)  in  a  slightly  different 
form,  namely  e  =  Cp^f-  This  different  form  is  due  to  the  use  of  a  bulk  drag  coefficient 

10  By  self-similarity  here  we  mean  that  K(z)  is  globally  dependent  on  U(z),  and  not  varying  in  its 
dependency  in  different  areas  of  the  canopy.  This  is  a  strong  assumption,  and  necessarily  approximate, 
but  facilitates  the  simple  but  powerful  analytical  model  developed  in  this  section. 
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Cd,  rather  than  a  series  of  coefficients  in  horizontal  slices  (Cd  varying  with  z).  By  setting 
X  =  U2(z )  and  (3  =  2e/ a,  the  differential  equation  reduces  to  the  simple  form 


_n 

dz2  H2X  ' 


(40) 


The  solution  of  this  equation  is  straightforward: 


X(z)  =  ci  sinh  +  c2  cosh  ^ 


H 


This  is  solved  in  conjunction  with  the  no-slip  condition  17(0)  =  0,  which  leads  to  c2  =  0. 
The  resulting  solution  for  U  becomes 


U(z)  =  uH 


/  sinh(/3z/i7) 
V  sinh(/3) 


z  <  H. 


(41) 


We  now  turn  our  attention  to  the  flow  profile  above  the  canopy,  ignoring  the  local 
scale  structural  complexity  of  the  roughness  sub-layer  (in  the  spirit  of  the  distributed  drag 
approach) ,  thereby  automatically  describing  this  layer  and  the  constant  stress  layer  above 
in  one  uniform  framework.  Following  [25],  we  generalise  (1)  for  velocity  above  the  canopy 
by  introducing  an  effective  canopy  displacement  height  d  and  roughness  thickness  zq- 


U(z)  =  —  \og(- — -V  z>H, 

7«  \  z0  J 

K(z)  =  7  Ku*(z  —  d),  z>H.  (42) 


The  diffusivity  profile  is  chosen  such  that  the  shear  stress  (35)  above  the  canopy  is  a  con¬ 
stant.  The  dimensionless  constant  7  is  a  deformation  parameter,  capturing  the  departure 
of  the  vertical  profiles  to  that  of  a  surface  without  a  canopy.  It  is  to  be  determined  either 
experimentally  or  on  physical  grounds  (see  [25,  26,  27];  7  ~  1.5  in  [25]).  We  expect  that  all 
our  equations  should  be  valid  (at  least  approximately)  for  a  value  7  =  1,  which  is  assumed 
below  unless  otherwise  specified.  This  is  in  keeping  with  our  philosophy  of  not  introducing 
extra  empirical  and  tunable  parameters  that  are  difficult  to  specify  in  a  robust  manner. 
However  it  is  clear  that  in  reality  gamma  will  vary  with  different  canopies,  and  this  will 
be  explored  in  future  work. 


In  a  similar  way  to  [25],  to  determine  the  unknowns  <7,  uh,  d  and  zq,  we  can  introduce 
four  physically  reasonable  boundary  conditions,  which  match  the  canopy  and  above-canopy 
profiles: 


U+{H)  =  U-(H), 


d  U+(z) 


dz 


dU-(z) 


z=H 


dz 


z=H 


K+(H)  =  K.(H), 


dK+(z ) 


dz 


d  K-(z) 


z=H 


dz 


z=H 


(43) 


Here  the  +  and  —  subscripts  denote  the  “upper”  and  “lower”  formulae,  valid  above  and 
below  z  =  H  respectively.  These  conditions  enforce  the  continuity  and  smoothness  of 
U{z)  and  K(z)  at  height  H.  The  last  condition,  imposing  smoothness  of  K(z)  is  new,  and 
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is  required  as  the  advection  diffusion  equation  (21)  used  to  find  the  plume  concentration 
explicitly  involves  the  first  derivative  of  K(z).11 

Solutions  (41)  and  (42)  are  matched  with  the  boundary  conditions  (43)  to  give  the 
following  relationships: 


d 

o 

H 

(7,t)2’ 

zo 

a 

H 

(7 K)2el  ’ 

UH 

1 

7«  ’ 

1%) 

1  far 

°) 

(7/c)2  V  2 

tanh  I  \l  —  = 

vv  aJ 

Equations  (44)  and  (45)  lead  to  the  useful  relationship 

d  +  zo 


H 


=  1  - 


a 

(7«) 


(1-0. 


(44) 

(45) 

(46) 

(47) 

(48) 


In  this  analysis  e  — »  oo  corresponds  to  an  infinitely  dense  canopy  and  e  — >  0  to  an  infinitely 
sparse  canopy. 

The  resulting  behaviour  of  the  velocity  profile  is  given  in  Figure  2  for  the  cases  of  a 
relatively  sparse  canopy  and  for  a  dense  canopy.  In  the  sparse  canopy,  the  wind  momentum 
penetrates  the  canopy  deeply,  so  the  drag  is  applied  directly  to  the  underlying  surface.  In 
the  opposite  case  of  a  dense  canopy,  very  little  wind  reaches  the  surface,  so  most  of  the 
drag  is  absorbed  by  the  canopy  elements.12  In  general  our  velocity  profile  is  similar  to  the 
one  found  by  DNS  in  [14]. 

It  is  readily  seen  that  as  e  — >  oo,  (47)  implies 


and  thus 


2(tk)2 


d  +  z0  —>  H, 
d  -»•  H, 
z0  —>  0. 


It  should  be  noted  that  this  limit  could  be  expected  to  be  identical  to  that  with  no 
canopy,  but  raised  by  height  H.  Substitution  shows  this  is  not  the  case,  and  is  an  artifact 
of  holding  the  velocity  constant  and  non  zero  at  the  canopy  top  (height  H)  while  matching 

11  It  should  be  noted  that  the  analysis  prior  to  the  imposing  of  these  boundary  conditions  is  identical  to 
that  by  Cowan  in  [38] .  However  Cowan  does  not  impose  smoothness  at  the  boundary.  It  is  the  addition  of 
the  expanded  boundary  conditions  of  the  first  derivatives  used  here  that  allows  the  determination  of  the 
key  flow  parameters  d  and  zq  in  terms  of  canopy  density  necessary  to  make  the  model  more  analytic. 

It  was  only  realised  after  the  fact  that  some  of  this  model  had  been  previously  derived  in  [38]. 

12This  is  validated  for  our  model  by  studying  the  shear  stress  at  the  ground  as  a  function  of  parameter 
e  (see  Figure  3). 
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Z 


H 


Figure  2:  Average  horizontal  velocity  profile.  The  pink/solid  line  is  for  a  sparse  canopy 
(e  =  0.5,  typical  urban)  and  the  blue/dashed  line  for  a  dense  canopy  (e  =  2,  forest). 


parameters  between  the  canopy  wind  profile  and  the  wind  profile  above  it.  This  shall  be 
explored  and  corrected  in  a  future  paper. 

In  the  opposite  limit  e  — >  0 

cr  ->  2(7k)2, 

d  +  z0  —>  2(e_1 -  1  )H, 
d  ->  -H, 
zq  — >  2He~1. 

The  main  feature  of  the  limit  e  — >  0  is  independency  of  the  parameters  d  and  zq  from 
e,  so  their  behaviour  is  universal.  Unfortunately,  the  skill  of  our  model  in  predicting  the 
limits  d(e  =  0)  and  zo(e  =  0)  are  not  very  advanced  and  needs  further  improvement  in 
these  extremes  (future  work  will  be  based  around  the  approach  of  Appendix  A.l).  It 
should  be  noted  that  the  expected  values  for  d(e  =  0)  and  zq(c  =  0)  are  0  and  the  viscous 
length  scale13,  respectively. 

The  plots  of  equations  (44)  and  (45)  showing  d  and  zq  as  functions  of  e  are  presented 
in  Figure  3a.  It  is  interesting  to  compare  these  plots  with  the  empirical  plots  and  formulas 
used  in  the  UDM — see  (11),  (12)  and  Figure  25  of  [1].  Since  e  is  proportional  to  A / 
the  agreement  between  the  results  for  e  >  1  (for  dense  canopies)  is  satisfactory.  In  the 
limit  of  the  sparse  canopy  e  <C  1,  the  agreement  is  not  visually  good,  since  this  limit 
involves  some  ambiguity  (for  a  discussion  see  Appendix  A.l).  However  away  from  the 
sparse  canopy  extreme,  the  U(z),  d  and  zq  profiles  are  quite  reasonable  and  agree  well 
with  experiment  (see  Section  4.4).  The  lower  limit  of  e  for  which  we  can  use  this  theory 
without  much  concern  is  determined  by  when  d  =  0.  The  corresponding  threshold  value 
of  e  is  e  ~  0.3,  which  is  smaller  than  typical  values  found  for  canopies  considered  in  this 
report  (see  below). 

13This  is  a  lower  bound  placed  upon  the  roughness  length,  as  viscosity  creates  effective  roughness  even 
for  perfectly  smooth  surfaces.  Usually  this  effect  is  swamped  by  the  presence  of  canopy  objects. 
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Relative  Heights 


Figure  3:  Left  -  Pink:  Relative  displacement  height  d/H  as  a  function  of  canopy  drag 
parameter  e  (solid  pink  line).  That  d  <  0  for  e  — >  0  is  an  artifact  of  the  functional  form 
imposing  a  behaviour  change  in  flow  at  height  H ,  even  for  an  infinitely  sparse  canopy. 
This  is  remedied  by  letting  H  — >  0  as  e  — >  0.  The  dashed  pink  line  is  a  sketch  of  d 
corresponding  to  the  UDM  (11).  Left  -  Blue:  Relative  roughness  height  zq  as  a  function 
of  canopy  drag  parameter  e.  Similarly  to  d,  the  fact  that  zo  f  0  for  e  — >  0  is  an  artifact  of 
the  functional  form  imposing  a  behaviour  change  in  flow  at  height  H ,  even  for  an  infinitely 
sparse  canopy.  This  is  remedied  by  letting  H  — ■>  0  as  e  — >  0.  The  dashed  blue  line  is  a 
sketch  of  zo  corresponding  to  UDM  (12).  Right:  Plot  of  wind  shear  stress  at  ground  level 
computed  with  the  proposed  model.  We  see  a  clear  monotonic  decrease  of  surface  shear 
stress  with  the  canopy  density  e,  as  intuitively  reasoned  on  physical  grounds. 


For  model  consistency  it  is  enlightening  to  study  the  shear  stress  at  ground  level,  given 
by  (35)  and  plotted  against  canopy  density  e  in  Figure  3b.  As  seen  in  this  plot,  for  very 
sparse  canopies,  the  shear  stress  is  concentrated  at  the  ground,  whereas  for  dense  canopies, 
the  wind  does  not  penetrate  to  the  ground  and  there  is  practically  no  shear  stress  there, 
as  is  required  on  physical  grounds.  Examining  Figure  3,  we  can  claim  that  for  our  model 
all  canopies  with  e  >  0.8  can  be  considered  dense. 

It  is  interesting  to  calculate  the  mixing  length  lm  used  by  other  models — e.g.  (18)  and 
(19)  discussed  in  Section  2.1 — based  on  our  framework.  The  formula  for  lm  is 


lm.  —  n 


U{z) 

U'(z) 


and  lm  is  plotted  for  sparse  and  dense  canopies  in  Figure  4.  This  can  be  compared  to  plots 
for  lm  found  in  the  work  of  Coceal  et  al  [14]  (Figure  21  of  that  paper).  The  general  trend 
of  our  profile  of  lm  matches  the  results  of  the  other  models  in  key  features,  such  as  that 
lm  increases  sharply  for  z  >  H,  and  is  non-zero  within  the  canopy. 


3.2  A  New  Model  for  the  Mean  Concentration  Profile 

The  new  model  presented  in  the  previous  section  for  flow  in  and  above  a  canopy,  combined 
with  the  universal  solutions  for  the  mean  concentration  profile  in  a  shear  flow  (24)  and  (30) 
lay  the  ground-work  for  a  new  model  of  concentration  profile  in  the  ABL  with  a  canopy. 
The  key  points  to  this  model  are: 
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Figure  4:  Plot  of  mixing  length  lm  for  our  flow  model.  The  pink  (long  dashed)  and  blue 
lines  are  for  e  =  1  and  4  respectively.  The  red  (short  dashed)  line  is  to  give  a  comparison 
of  the  gradient  of  Figure  21  in  [14]  for  z  >  H. 


1.  Since  the  mean  flow  above  the  canopy  (42)  can  be  modelled  with  the  two  parameters 
d  and  zo,  expressed  in  Section  3.1  in  terms  of  canopy  density  e  and  height  H  through  (44), 
(45)  and  (47),  the  power-law  velocity  profile  (8)  with  the  same  value  of  these  parameters 
is  also  a  realistic  model  of  the  mean  flow  in  the  canopy14.  For  this  power  law  profile, 
the  stretched  exponential  solutions  for  concentration  (24)  and  (30)  are  exact  solutions 
of  the  turbulent  diffusion  equation  (21),  so  we  come  to  the  conclusion  that  the  formula 
(30)  should  be  a  realistic  model  of  concentration  distribution  for  z  >  d,  if  (  is  replaced 
by  C  =  (z  ~  d)/zo,  where  d  is  the  height  that  the  canopy  has  displaced  the  plume  by, 
and  is  defined  by  our  flow  model  (see  Figure  1).  We  acknowledge  that  the  new  model  for 
canopy  flow  derived  in  the  previous  section,  which  uses  a  log  law  profile,  is  mathematically 
inequivalent  to  the  power  law  profile  employed  to  derive  the  vertical  mean  concentration 
distribution,  and  thus  the  models  are  in  some  sense  incompatible.  However  as  shown  in 
Section  2.1,  within  the  ABL  the  two  models  are  similar  enough  to  be  interchangeable. 
Thus  we  take  a  practical  approach  here  that  employing  these  different  velocity  functions 
for  different  purposes  is  the  only  way  to  make  progress  in  deriving  simple  analytical  models 
in  the  first  instance. 

What  is  the  concentration  distribution  for  z  <  d?  The  simplest  assumption  is  to 
assume  that  far  from  the  source  the  flow  in  the  canopy  is  well-mixed  and  is  in  dynamical 
equilibrium  (all  concentration  gradients  have  been  smoothed  out  by  turbulent  diffusion). 
This  leads  to  the  conclusion  that 

Cz(x,  z)  «  C(x,d),  z  <  d.  (49) 

This  “clipped  stretched  exponential”  profile  is  a  good  approximation  everywhere  outside 

14There  is  a  cost  to  taking  this  equivalence,  in  that  the  power  law  profile  is  not  universal,  and  introduces 
constants  which  will  need  to  be  fitted  by  comparison  to  experiment,  however  this  a  traded  off  against  a 
functional  form  of  the  velocity  profile  for  which  analytic  plume  concentration  calculations  are  possible. 
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the  the  so-called  source  region,  which  is  typically  one  to  two  obstacle  rows  downstream 
in  the  obstacle  arrays  considered  in  this  report  and  other  similar  studies.  A  supporting 
analytical  argument  can  be  found  in  Appendix  A. 2. 

2.  An  alternative  model  to  the  “clipped  stretch  exponential”  model  (49)  can  be  for¬ 
mulated  in  terms  of  the  power- law  exponent  a  in  (24). 

In  Appendix  A. 2  it  is  shown  that  near  the  bottom  of  the  canopy  the  advection-diffusion 
equation  can  be  solved  to  give  a  mean  concentration  profile  with  a  =  2.  Thus  a  smooth, 
continuous  function  for  the  concentration  profile  can  be  arrived  at  by  assuming: 

a  =  ao(l  —  Aexp(— z/d)),  (50) 

A  =  1  —  2/ao 

and  inserting  this  into  mean  concentration  profile  (24).  The  ansatz  for  a  has  been  designed 
to  change  behaviour  at  the  effective  canopy  height  d. 

Conceptually,  close  to  the  ground  this  model  predicts  that  the  effect  of  the  canopy  is 
to  restrict  the  wind  and  to  cause  the  tracer  to  disperse  in  a  Gaussian  manner.  Then  above 
the  canopy,  the  dispersal  is  dictated  by  the  power  law  behaviour  of  the  wind  and  diffusivity 
functions.  In  Section  4.5  this  model  is  validated  against  our  experimental  data-set  and 
provides  a  close  match. 

It  should  be  emphasised  that  the  vertical  dependency  of  a  introduced  in  (50)  will 
modify  not  only  the  power-law  exponent  in  concentration  formulae  (24)  and  (30),  but  also 
the  pre-exponential  term  Co,  which  contains  a  within  xq.  From  physical  grounds,  such 
complex  modification  of  the  concentration  profile  is  imposed  by  the  self-similarity  property 
of  the  solutions  (24)  and  (30). 

3.  The  next  new  feature  of  our  model  is  an  “effective  elevation”  of  the  tracer  source, 
i.e.  non-zero  value  of  the  parameter  £0  in  (30),  even  for  ground- level  sources.  The  physical 
reasoning  behind  this  is  an  observed  effect  of  the  mean  flow  being  pushed  up  by  the 
canopy.  This  is  coupled  with  the  effect  of  the  highly  turbulent  recirculation  region  behind 
the  obstacle  where  the  source  is  placed,  rapidly  mixing  the  tracer  to  fill  the  entire  wake  up 
to  the  canopy  height  with  contaminant.  This  effect  is  most  pronounced  when  the  source  is 
placed  just  in  front  of  a  canopy  object.  Such  effects  have  been  observed  in  many  previous 
obstacle  dispersion  experiments  [7,  12,  28,  34,  36].  The  experimental  evidence  for  this 
with  the  currently  considered  data-set  is  covered  in  Section  4.5. 

The  resulting  pattern  of  the  pollutant  is  as  if  it  were  released  from  an  elevated  source 
(Figure  13  below  gives  an  example  of  this  effect  from  our  data).  Including  a  non-zero 
parameter  £0  in  our  model  for  ground-level  sources  permits  a  realistic  description  of  the 
vertical  variability  of  the  concentration  profile  in  the  “near  source”  field  before  the  con¬ 
centration  distribution  approaches  the  “clipped  stretched  exponential”  regime  described 
above. 

It  is  possible  to  provide  a  simple  estimate  of  the  value  of  Co  in  terms  of  the  relative 
source  position  and  the  averaged  canopy  height  H  for  ground  level  sources  within  the 
canopy.  It  is  evident  that  0  <  Co  5-  G  s°  as  the  first  approximation  we  can  take 

Co  =  (77  -  az0)/H,  (51) 
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Figure  5:  Illustration  of  the  concept  of  a  virtual  source  induced  by  the  wake  regions  of 
the  obstacle  array.  The  virtual  source  is  lifted  to  height  £o H  (52). 


for  Co  <  1  and  Co  =  0  otherwise  (see  Figure  5).  The  plume  spread  aZQ  is  the  corresponding 
vertical  dispersion  parameter  taken  if  no  canopy  were  present,15  and  can  be  estimated 
from  the  well-known  kinematic  condition  [4]  az o  ~  V±x/V\\  ~  u*x/Uq  ~  Cpx,  where  rt*  is 
the  friction  velocity,  Uq  is  the  undisturbed  velocity  far  from  the  surface  (at  the  top  of  the 
boundary  layer)  and  Cf  is  the  drag  coefficient  for  the  flow  over  a  flat  smooth  surface  [4] . 
Thus  we  obtain 

Co  =  (H-  ipCfXH)/H,  (52) 

where  Xu  is  a  distance  from  the  source  to  the  first  canopy  element  downstream  and  if  is 
a  numerical  factor  to  be  determined. 

To  summarise  the  urban  plume  model  discussed  above,  we  have  Cy(x,y )  unchanged 
from  (23),  and  combining  the  elements  above,  we  obtain  the  following  model  for  Cz(x,  z ): 


Cz 

Cb(s,C) 

C(*) 

Co 


(  C0(x,z)exp(-B(x)(Co  +  Cq)K-i/(25(x)(CCo)q/2),  \iz>d 

|  Cz(d),  if  0  <  z  <  d, 

— ,  B(x)  = 

z0  X 

H  -  ifCfXH  _  uqZq  (1  +  to) 

II  azI\o  a 


(53) 

(54) 

(55) 

(56) 


Including  the  varying  a  hypothesis  instead  of  a  having  constant  behaviour  below  d  is  also 
possible,  but  leads  to  complications  that  will  be  resolved  in  future.  Instead,  in  the  next 
section,  varying  a  is  studied  without  an  effective  source  rise. 

The  verification  of  the  components  of  this  model  is  discussed  in  the  following  section. 

15Thus  not  to  be  confused  with  the  actual  plume  spread  derivable  from  the  canopy  model  being  discussed. 
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4  Data  Analysis 


4.1  Introduction 

The  plume  dispersion  models  described  in  this  report  do  not  give  self-contained  predictions 
of  the  contaminant  concentration  fields.  Instead  they  predict  the  functional  form  of  the 
concentration  distributions  based  on  assumptions  about  the  nature  of  the  underlying  flow. 
Analysis  of  experimental  data  can  validate  the  assumptions  about  the  flow  used  to  derive 
the  dispersion  models.  The  data  can  also  fix  the  free  parameters  of  the  dispersion  models 
to  give  predictive  power  to  the  models. 

The  data  presented  here  were  collected  in  water  channel  simulations  carried  out  in  the 
Coanda  Research  &  Development  Corporation’s  (Burnaby,  British  Columbia,  Canada) 
environmental  water  channel.  The  simulations  were  conducted  in  collaboration  with  De¬ 
fence  Research  and  Development  Canada  (DRDC),  which  had  a  multi-year  program  in 
the  development  of  the  channel  and  collection  of  extensive  velocity  and  1-dinrensional 
(ID)  tracer  dispersion  data-sets.  DSTO  extended  the  experimental  program  to  investi¬ 
gate  2-dimensional  (2D)  sampling  of  the  tracer  concentration,  as  well  as  in-canopy  velocity 
measurements,  as  described  below.  The  water  channel  is  specially  designed  for  dispersion 
modelling,  and  has  a  working  section  of  10  m  length  with  a  cross  section  of  1.5  m  width 
and  0.9  m  height.  We  will  give  a  brief  description  of  the  water  channel  experiment  here, 
to  establish  the  context  behind  the  data.  Considerably  greater  details  of  the  water  chan¬ 
nel  experiments  and  the  instrumentation  deployed  are  available  in  a  report  provided  by 
Coanda  [29],  and  previously  published  work  using  this  data  [28,  30]. 


4.2  Brief  Experimental  Description 

The  upwind  part  of  the  working  section  of  the  water  channel  was  used  to  generate  an 
appropriately  scaled  model  of  the  neutrally-stratified  ABL.  A  thick  rough- wall  boundary 
layer  in  the  water  channel  of  about  0.3  nr  thickness  was  generated.  The  combination  of  saw¬ 
tooth  fence  and  turbulence  grid  made  of  square  bars  19x19  mm  was  used  to  immediately 
produce  some  of  the  larger-scale  turbulent  eddies  in  the  boundary  layer.  Downstream  of 
the  saw-tooth  fence,  the  floor  of  the  water  channel  was  covered  with  a  black  anodised 
expanded  metal  mesh  of  height  4  mm  with  a  total  length  of  6  nr  to  give  a  long  fetch  of 
surface  roughness  for  the  upstream  flow  development.  The  flow  then  encountered  either 
further  nresh  for  the  case  of  open  terrain  with  no  obstacles,  or  the  model  obstacle  array, 
which  was  mounted  on  a  turntable  that  extended  1.25  nr  in  the  streanr-wise  direction. 
Downstream  of  the  obstacle  array,  the  flow  encountered  a  section  covered  with  the  same 
black  anodised  expanded  nretal  mesh  that  was  used  for  the  initial  upstream  approach  flow 
development. 

Two  approaches  were  used  to  measure  velocity  statistics  within  the  boundary  layer 
of  the  water  channel.  A  dual-beam  laser  doppler  velocinreter  (LDV)  fibre-optic  probe 
was  used  in  the  original  experiments  commissioned  by  DRDC.  To  complement  this  data 
set,  DSTO  commissioned  some  acoustic  doppler  velocinretry  (ADV)  measurements  using 
a  Sontek  MicroADV  instrument.  Both  these  instruments  have  their  strengths  and  weak¬ 
nesses  [31],  and  a  combination  of  the  two  data  sets  provides  a  complete  and  complementary 
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set  of  velocity  measurements.  Vertical  profiles  of  mean  wind  and  other  turbulence  param¬ 
eters  were  measured  in  the  equilibrium  (fully-developed)  boundary  layer  upstream  of  the 
obstacle  arrays.  The  height  of  the  boundary  layer,  taken  to  be  the  height  where  the  mean 
wind  was  99%  of  the  free-stream  value,  was  found  to  be  275  mm.  At  this  point,  the  mean 
wind  speed  was  0.375  nr  s-1. 

For  the  water  channel  simulations,  a  ground-level  point  source  was  used,  emitting 
sodium  fluorescein  dye  tracer  at  a  constant  flow  rate.  An  example  of  one  particular  obstacle 
array  is  presented  in  Figure  6  (left  panel).  The  source  was  located  at  the  first  street  canyon, 
and  the  measurements  were  conducted  for  the  various  models  of  urban  terrain  (obstacle 
array  configurations).  The  size  of  the  “atomic”  cubic  element  of  the  various  obstacle  arrays 
was  77  =  31.75  mm.  A  number  of  Urban  Arrays  of  different  heights  (1 77,  2 77,  3 77  and 
their  random  combination)  and  different  space  patterns  (regular  and  random)  were  used  in 
the  simulations.  Measurements  in  the  so-called  “Tombstone  Array”  were  also  conducted 
to  enable  comparison  to  a  similar  experiment  conducted  using  ID  line-scans  by  DRDC. 
The  Tombstone  Array  was  constructed  from  vertical  stainless  steel  tabs  (60x10x1  mm) 
installed  vertically  into  laser  cut  slots  and  arranged  in  staggered  formation.  For  Urban 
Arrays,  measurements  were  conducted  at  six  downstream  locations  (2.077,  4.077,  1077, 
1677,  2677),  and  for  the  Tombstone  Array  at  five  downstream  locations  (Is,  2s,  4s,  7s, 
12s,  where  s  =  44  mm). 


AnayOOl 


Cell  at  Row  2 


0.00 

-7.94 

-15.88 

-23.81 

C 
O 

■|  -31.75 

CL 
> 

-39.69 

-47.63 

-55.56 

-63.50 

317.50  325.44  333.38  341.31  349.25  357.19  365.13  373.06  381.00 

X  position 

Figure  6:  Left:  Coanda  Urban  Array  001  consisting  of  regularly  spaced  cubic  objects 
with  packing  density  0.25.  Right:  Location  of  positions  within  a  cell  unit  where  velocity 
measurements  were  taken.  The  dark  shaded  area  indicates  the  position  of  the  object  within 
the  unit  cell. 


The  instantaneous  concentration  field  in  the  dispersing  dye  plume  was  measured  using 
the  2D  laser-induced  fluorescence  (LIF)  technique,  which  permitted  simultaneous  multi¬ 
point  concentration  measurements  to  be  acquired  with  high  spatial  and  temporal  reso¬ 
lution.  A  fluorescent  dye  is  emitted  from  a  source,  and  then  disperses  in  the  turbulent 
boundary  layer.  The  dye  encounters  a  laser  sheet  that  causes  the  dye  to  fluoresce.  Each 
pixel  of  the  real  time  fluorescing  image  can  be  accurately  calibrated  to  a  concentration 
field  and  captured  in  a  CCD  camera.  Each  2D  scan  consisted  of  696x260  pixels,  giving  a 
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spatial  resolution  of  about  1.4  mm.  The  2D  scans  were  sampled  for  1000  s  at  23  frames 
per  second,  giving  data  sets  of  about  23,000  images  per  experiment. 

In  addition  to  low  level  data,  Coanda  has  also  supplied  analysis  programs  and  anal¬ 
ysed  data.  A  software  framework  has  been  developed  here  at  the  DSTO  to  allow  the 
comparison  of  plume  dispersion  models  and  the  Coanda  data.  This  software  framework 
utilises  the  lower-level  Coanda  software  and  processed  data,  and  applies  higher-level  anal¬ 
ysis.  For  preliminary  data  analysis  the  software  was  developed  in  MATLAB.  A  more 
advanced  framework  is  written  in  the  C+- 1-  programming  language.  It  provides  a  conve¬ 
nient  environment  where  large  amounts  of  processed  data  can  be  collected  and  analysed. 
Software  to  read  and  analyse  low  level  Coanda  data  has  also  been  developed.  Fitting  and 
other  analyses  are  performed  using  the  ROOT  analysis  package  [32]  developed  by  the  high 
energy  physics  community. 


4.3  Coanda  Dataset — Preliminary  Analysis 

The  Coanda  data-set  has  been  used  in  several  studies  of  plume  dispersion  [28,  30,  33]. 
In  this  report  the  2D  Coanda  data-set  is  used  to  validate  analytical  predictions  of  the 
models  developed  in  the  previous  section.  As  a  first  step  for  such  validation,  we  verify  the 
existence  of  universal  scaling  (29)  and  existence  of  stretched  exponential  profiles  (24)  and 
(30)  in  a  complex  canopy. 

According  to  (29),  the  ratios  Oy/n  and  az/ /x  should  not  depend  on  the  downstream 
position.  This  fact  was  demonstrated  in  a  previous  paper  [30]. 

Also,  that  the  cross-stream  mean  concentration  is  given  by  Gaussian  distribution  (23) 
has  been  observed  in  many  previous  studies  [12,  28,  33,  34,  36].  This  Gaussian  behaviour 
can  be  seen  in  figures  11  and  13. 

The  stretched  exponential  solution  (24)  was  used  to  match  the  mean  concentration 
profile  in  the  vertical  direction  for  the  turbulent  surface  layer  above  both  smooth  and 
rough  surfaces.  The  results  are  presented  in  Figure  7.  As  can  be  seen  the  stretched 
exponential  model  is  consistent  with  the  measured  concentration  data. 

In  Figure  8  the  plot  of  y  =  \og{\og(Cmax/C))  against  x  =  log  z  is  presented  (smooth 
surface  case) ,  so  each  stretched  exponential  profile  corresponds  to  a  straight  line  (with  the 
slope  of  such  a  line  determined  by  the  power  exponent  a  in  (24)).  From  the  left-hand 
plot,  we  observe  that  the  stretched  exponential  function  is  a  plausible  model  for  the  mean 
concentration  profile  in  the  vertical  direction  for  the  case  of  a  simple  sheared  boundary 
layer  (all  plots  collapse  under  normalised  coordinates). 

The  righthand  plot  in  Figure  8  shows  the  effect  of  the  regular  cubical  obstacle  array 
on  the  standard  vertical  profile  form.  We  can  see  that  the  canopy  causes  a  deformation  of 
the  slope  around  height  H,  i.e.  plot  of  log(log(C)  becomes  horizontal  within  the  canopy. 
Because  of  mixing  and  flow  homogenisation,  the  value  of  a  tends  to  0  in  this  region.  This 
is  consistent  with  the  “clipped  stretched  exponential  solutions”  (49)  and  (53).  However 
we  do  not  use  a  — >  0  to  model  the  concentration  profile  using  a  varying  a  formalism  (as 
in  (50))  as  this  contradicts  the  result  of  Appendix  (A. 2). 

From  these  plots  we  estimate  a  value  of  the  a  exponent  of  the  stretched  exponential 
distribution  for  various  downstream  positions.  We  find  a  good  agreement  between  values 
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Figure  7;  Vertical  cross-sections  of  the  normalised  mean  concentration  as  a  function  of 
downwind  position.  Solid  line  is  the  stretched  exponential  fit  (24).  Ca  depicts  the  time 
averaged  concentration  at  that  height  and  downstream  position,  Ca,max  is  the  maximum 
concentration  along  this  vertical  axis.  aaz  is  the  standard  deviation  in  the  vertical  position. 
Description  of  the  obstacle  arrays  can  be  found  in  [29]. 
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Figure  8:  Plot  of  log(log(Cmax/Ca))  against  normalised  vertical  coordinate  log z/aaz. 
The  slopes  of  curves  are  determined  by  the  power  exponent  a  in  (2f). 


of  a  found  from  data  presented  in  Figure  8.  (and  plotted  in  Figure  9)  with  our  theoretical 
estimates  (ac  =  1  +  2m  =  1.54,  a«o  =  1  +  m  =  1.27)  for  a  measured  value  of  m  =  0.27 
(obtained  from  fitting  to  data  in  Figure  11). 


4.4  Coanda  Flow  Velocity  Measurements 

The  velocity  of  the  atmospheric  flow  within  and  above  an  urban  canopy  is  the  key  driver  of 
plume  dispersion.  For  the  models  described  in  this  report  the  mean  stream-wise  velocity 
and  its  height  dependence  are  important  parameters. 

The  Coanda  data-set  includes  measurements  of  the  velocity  fields  over  a  range  of 
mock  urban  canopies.  Individual  data  points  give  a  measurement  of  two  of  the  velocity 
components  at  a  single  point  in  the  flow.  With  repeated  measurements,  the  full  3-D  flow 
structure  can  be  characterised  at  a  range  of  heights  and  positions  within  and  above  the 
canopy. 

Figure  10  shows  an  example  of  the  velocity  information  that  can  be  extracted  from 
the  Coanda  data-set.  The  configuration  of  the  obstacle  array  for  this  measurement  is 
designated  Urban  Array  001.  This  array  consists  of  a  regular  pattern  of  cubic  objects  as 
shown  in  Figure  6.  Also  shown  in  Figure  6  (right  panel)  are  the  positions  within  one  of  the 
array  sub-units  where  the  velocity  measurements  of  Figure  10  were  taken.  The  velocity 
measurements  in  Figure  10  show  that  the  flow  velocity  above  the  canopy  is  uniform  over 
the  array  (ie  the  data  points  all  follow  a  common  line).  That  the  data  points  within  the 
canopy  have  a  significant  spread  in  velocities  indicates  that  there  are  complex  flow  patterns 
depending  on  the  location  of  the  measurement  with  respect  to  the  objects  of  the  array.16 

16It  should  be  noted  that  the  measurement  points  are  not  ideally  placed  for  this  fit  to  give  a  true  average 
of  the  velocity  profile  of  the  system.  Rather  this  plot  gives  a  representative  idea  of  the  behaviour  of  the 
velocity  profile.  Negative  velocities  correspond  to  eddies  directly  behind  the  array  objects. 
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Figure  9:  Variation  in  the  parameter  a  of  (24)  as  a  function  of  downstream  distance 
from  the  source  for  the  Coanda  data  with  no  obstacles. 


Figure  10:  Measurements  of  height  z  vs.  mean  stream-wise  velocity  U  for  the  Coanda 
Urban  Array  001.  The  data  points  come  from  measurements  at  different  representative 
points  around  a  single  array  object.  The  configuration  of  the  array  and  the  location  of 
measurement  points  can  be  seen  in  Figure  6  (right  panel).  The  dashed  blue  line  indicates 
the  height  of  the  objects  in  the  array.  The  solid  line  is  a  fit  of  all  of  the  data  to  the 
functional  form  given  in  (41)  and  (42). 
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4.5  Coanda  Mean  Concentration  Measurements 

As  mentioned  in  the  Introduction,  the  long  term  goal  of  the  work  presented  in  this  report  is 
to  find  models  that  can  describe  the  concentration  fluctuations  in  a  plume.  An  important 
first  step  in  developing  these  models  is  to  describe  the  time  averaged  behaviour  of  the 
plume.  The  basic  form  of  the  stretched  exponential  solution  for  mean  concentration  in  a 
plume  is  shown  in  (24). 

Underlying  this  solution  is  the  description  of  the  flow,  namely  that  the  mean  down¬ 
stream  velocity,  U,  and  diffusivity,  K  vary  with  with  height  as  a  power  law.  Equation  (24) 
has  four  free  parameters  associated  with  the  flow:  U  and  K  at  a  reference  height  and 
their  respective  power  law  height  indexes  m  and  n.  If  the  absolute  concentration  is  not 
normalised  this  gives  a  total  of  five  free  parameters.  Clearly  the  number  of  free  parameters 
is  a  major  difficulty  in  validating  (24).  Part  of  the  data  analysis  work  described  here  has 
been  to  maximise  the  use  of  available  data  to  constrain  the  fit  parameters. 

For  the  Coanda  data-set,  the  measured  concentrations  are  normalised  to  a  known 
source  strength.  Using  this  information  it  is  possible  to  fix  the  absolute  concentration 
normalisation  of  (24)  rather  than  have  it  as  a  fitted  parameter.  Further  constraints  can 
be  put  on  fitted  parameters  by  demanding  consistency  in  the  parameters  of  the  flow  at 
different  locations.  This  has  been  achieved  by  performing  a  global  fit  for  flow  parameters 
over  plume  concentration  profiles  at  different  downstream  distances  from  the  source. 

Figure  11  shows  an  example  of  fitting  (24)  to  Coanda  data  for  the  case  where  there 
is  no  urban  canopy.  The  figure  shows  that  if  the  fit  parameters  are  allowed  to  vary  for 
each  profile,  the  functional  form  of  (24)  can  provide  an  excellent  match  to  the  data.  If  a 
global  fit  is  done,  demanding  the  same  flow  parameters  in  all  parts  of  the  flow,  then  the 
fit,  although  not  perfect,  is  still  a  very  good  description  of  the  data. 

Figure  12  shows  fits  of  a  stretched  exponential  form  to  Coanda  data  where  an  urban 
array  is  present  (Urban  Array  001).  The  presence  of  the  canopy  causes  a  vertical  dis¬ 
placement  in  the  flow  as  discussed  in  Section  3.2.  In  Figure  12  the  vertical  concentration 
profiles  have  been  fit  to  the  form  given  in  (53),  with  displacement  height  d  calculated  from 
(44)  to  be  11  mm.  Within  the  canopy  the  concentration  is  assumed  to  be  constant.  If  the 
flow  parameters  are  allowed  to  vary  at  different  locations  in  the  array  then  the  functional 
form  of  (53)  describes  that  data  well. 

If  a  global  fit  is  done,  with  consistent  flow  parameters  throughout  the  array,  then 
the  fit  works  poorly.  The  reason  for  this  is  probably  that  while  the  functional  form  (53) 
does  predict  that  the  amount  of  material  in  a  plane  perpendicular  to  the  direction  of  flow 
attenuates  with  distance  from  the  source17,  the  concentration  attenuation  in  the  Coanda 
experiment  is  much  too  high  to  be  explained  by  this  effect  (the  functional  form  predicts  an 
8-10%  loss  between  the  508  mm  and  the  825  mm  points,  whereas  what  is  lost  is  closer  to 
50%).  Rather  material  is  apparently  being  lost  because  the  sensors  cannot  detect  a  local 
concentration  below  a  specific  threshold,  and  thus  when  a  plume  fluctuation  drives  the 
concentration  too  low,  especially  when  the  average  concentration  in  the  region  is  already 
low,  then  the  sensor  records  that  no  tracer  is  present.  Alternative  explanations  are  that 
the  flow  is  non  homogeneous  over  the  array  (quite  possible  as  the  array  designers  put 

17This  is  because  conservation  of  concentration  flux  at  the  higher  wind  speeds  of  greater  altitudes 
corresponds  to  lower  concentrations. 
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Figure  11:  Plume  concentration  profiles  from  Coanda  data:  a  ground  level  release  over 
flat  terrain  (no  canopy).  The  plots  on  the  left  show  the  vertical  concentration  profile  at 
different  downstream  distances  (black  and  blue  points).  Concentration  has  been  normalised 
by  concentration  of  dye  in  the  source  pipe  fluid  and  is  dimensionless.  The  plots  on  the 
right  show  the  horizontal  plume  profile.  In  each  plot  the  magenta  line  shows  the  results  of 
fitting  (24)  (vertical  profiles)  and  a  Gaussian  function  (horizontal  profiles)  while  allowing 
the  flow  parameters  to  vary  at  each  downstream  position.  The  red  and  green  lines  in  the 
vertical  profile  plots  show  the  form  of  (24)  if  the  fit  parameters  are  derived  from  a  global 
fit  to  all  vertical  concentration  profiles. 
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Figure  12:  Same  plots  as  shown  in  Figure  11,  but  for  plume  dispersion  over  Coanda 
Urban  Array  001.  Vertical  axis  is  concentration  and  has  been  normalised  by  concentration 
of  dye  in  the  source  pipe  fluid  and  is  dimensionless.  See  Figure  6  for  details  of  this  array 
configuration.  The  source  has  been  located  at  ground  level  in  a  canyon  between  rows  of 
objects.  The  functional  form  used  in  the  fits  is  given  by  (2j).  The  dashed  lines  in  the  left 
hand  plots  refer  to  a  global  fit  over  all  the  data. 
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the  plume  sources  near  the  beginnings  of  the  arrays,  so  turbulent  conditions  are  likely  to 
be  different  at  the  source  than  near  the  middle)  or  that  the  model  is  missing  some  key 
behaviour.  Further  investigation  is  required  to  resolve  this. 

The  vertical  displacement  of  the  flow  caused  by  a  canopy  is  an  example  of  a  global  effect 
of  a  canopy  on  the  flow.  In  the  Coanda  data  we  can  also  see  local  effects  on  the  plume 
due  to  structures  in  the  canopy.  Figures  12  and  13  show  concentration  profiles  taken 
under  identical  flow  and  array  configurations.  The  only  difference  is  that  in  Figure  12 
the  source  was  placed  in  a  “canyon”  between  rows  of  objects  in  the  canopy,  while  in 
Figure  13  the  source  was  placed  directly  upstream  of  an  object.  The  vertical  concentration 
profile  in  Figure  13  shows  a  pronounced  peak  that  is  higher  than  that  seen  in  Figure  12. 
We  hypothesise  that  the  local  structure  in  front  of  the  source  is  pushing  part  of  the 
plume  upward  simulating  a  non  ground-level  source.  This  type  of  effect  has  been  observed 
consistently  in  previous  obstacle  array  experiments  [12,  28,  34,  36].  To  fit  the  data  we  have 
used  the  form  of  the  raised  clipped  stretched  exponential  (53)  and  we  have  used  (44)  to 
calculate  the  effective  canopy  height  d=  11  mm  and  made  an  ad-hoc  assignment  of  source 
height  of  10  mm 18 .  As  can  be  seen,  the  fit  of  the  functional  form  to  the  experimental  data 
is  currently  faulty.  While  the  functional  form  can  indeed  produce  the  same  flat  region 
within  the  canopy  followed  by  a  rise  in  concentration  above  the  canopy,  which  fades  with 
distance,  the  algorithm  used  in  obtaining  this  plot  has  not  yet  been  adjusted  to  account  for 
the  effective  loss  of  tracer  as  discussed  for  Figure  12.  Initial  calculations  suggest  that  once 
this  weakness  in  the  experimental  data  is  accounted  for,  (53)  should  fit  the  experimental 
data  well. 

The  plot  in  Figure  14  represents  a  fit  with  the  varying  a  of  equation  (50).  We  observe 
that  this  model  provides  a  satisfactory  fit  in  the  region  close  to  the  source,  although  slightly 
worse  than  for  the  “clipped  stretched  exponential”  in  Figure  12.  As  expected,  the  two 
different  functional  forms  and  make  mostly  no  difference  above  the  canopy. 


5  Conclusions  and  Future  Work 


The  high  fidelity  modelling  of  tracer  dispersion  in  heterogeneous  canopies  is  a  compu¬ 
tationally  intensive  task  that  still  requires  long  execution  time  and  significant  computer 
resources.  In  the  current  report  we  have  shown  a  new  approach  for  development  of  simpli¬ 
fied  flow  and  dispersion  models  for  urban  environments  that  are  easily  treated  numerically 
in  an  operational  environment  with  minimal  computer  resources.  The  framework  estab¬ 
lished  does  not  intend  to  offer  a  new  generation  of  dispersion  models  to  replace  existing 
operational  models,  but  rather  provides  a  link  between  empirical,  fast  models  such  as  the 
UDM,  and  computationally  intensive  CFD  models.  The  models  can  thus  justify  some 
heuristic  approximations  of  the  UDM,  and  provide  an  important  performance  check  for 
CFD  models.  This  enables  a  more  comprehensive  understanding  of  the  problem,  with 
greater  clarity  in  the  approximations  and  assumptions  made  for  various  dispersion  mod¬ 
els.  We  have  presented  experimental  results  from  water  channel  experiments  to  support 
the  theoretical  findings. 

1SA  best  guess  assignment  is  used  at  this  stage  as  the  problems  with  the  current  global  fit  algorithm 
means  ip  in  (52)  is  yet  to  be  determined. 
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Figure  13:  Identical  plots  and  array  configuration  as  for  Figure  12  except  the  source  has 
been  placed  directly  behind  an  object  in  the  array.  The  functional  form  used  in  the  fit  is 
given  by  (53).  Vertical  axis  is  concentration  and  has  been  normalised  by  concentration  of 
dye  in  the  source  pipe  fluid  and  is  dimensionless. 
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Figure  1 4-  Example  of  model  fit  for  ground  source  in  a  canopy  array.  Vertical  axis  is 
concentration  and  has  been  normalised  by  concentration  of  dye  in  the  source  pipe  fluid  and 
is  dimensionless.  In  each  plot  the  magenta  line  shows  the  results  of  fitting  (24)  (vertical 
profiles)  and  a  Gaussian  function  (horizontal  profiles),  but  with  a  from  (50).  Better  fit  in 
the  source  region  (top  plot). 
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The  framework  is  also  important  as  a  first  step  in  the  development  of  models  of  con¬ 
centration  fluctuations  in  urban  environments.  Such  models  facilitate  the  generation  of 
realistic  plumes  for  synthetic  environments  used  in  operational  analysis  studies.  Thus  the 
strong  internal  variability  and  intermittency  of  the  concentration  field  within  a  dispersing 
plume  will  be  able  to  be  captured.  This  will  enable  the  calculation  of  realistic  fluctuat¬ 
ing  CBR  challenge  levels,  which  will  facilitate  studies  on  optimal  detection  systems  for 
CBR  threats,  provide  the  necessary  input  for  models  of  toxicity  response  and  data  fu¬ 
sion  algorithms  for  enhanced  situational  awareness,  and  is  a  cornerstone  for  the  rigorous 
quantification  of  uncertainty  in  dispersion  predictions. 

Scope  exists  to  refine  the  proposed  modelling  framework.  This  includes  resolution  of 
the  un-physical  behaviour  of  the  displacement  height  d  as  the  canopy  drag  parameter 
e  becomes  very  small.  To  provide  a  truly  unified  framework,  a  single  functional  form 
for  the  velocity  framework  would  also  be  preferable.  The  problem  here  is  that  there  is 
no  analytical  solution  for  the  advection  diffusion  equation  in  a  log-law  boundary  layer — 
the  traditionally  preferred  description  of  a  sheared  flow  boundary.  With  the  stronger 
theoretical  underpinning  of  the  power  law  velocity  profiles  due  to  Barenblatt  [5,  6],  this 
could  be  a  preferred  option,  however  the  flow  model  described  in  Section  3.1  would  be 
ill-determined  with  the  use  of  this  functional  form.  There  are  options  for  a  more  complex 
specification  of  the  profile  with  a  variable  power  law  exponent  a  that  would  enable  a 
solution  here,  but  this  added  complexity  and  extra  parameterisation  was  not  pursued  in 
the  first  instance  as  reported  here. 

The  equation  (53)  currently  models  what  should  happen  to  the  tracer  in  the  plume 
spreading  over  an  urban  canopy.  However  when  comparing  this  model  to  experimental 
data,  a  model  of  how  the  equipment  ‘sees’  a  particular  turbulent  plume  concentration 
is  also  needed,  in  order  to  accurately  compare  experiment  and  theory.  This  is  because 
any  experimental  setup  is  not  perfect.  In  the  case  of  the  Coanda  experiment,  a  weakness 
lies  in  the  equipment  not  being  able  to  sense  a  concentration  below  a  critical  threshold, 
and  as  such  part  of  the  plume  is  not  captured  by  the  experiment,19  resulting  in  apparent 
non-conservation  of  the  plume  in  the  data  and  deformation  of  the  plume  shape.  Such  a 
step  is  essential  for  the  model  to  be  able  to  predict  not  only  what  a  plume  would  do,  but 
how  sensors  and  experiment  would  see  the  plume.  It  would  also  allow  completion  of  the 
validation  of  the  model  with  the  Coanda  data-set,  such  as  in  Figures  12  and  13. 

This  is  suspected  to  be  a  key  reason  why  a  global  fit  of  parameters  for  the  model  over 
the  whole  scope  of  a  plume  affected  by  an  underlying  canopy  is  poor.  Other  possible 
reasons  are  that  the  COANDA  experiment’s  flow  may  be  non  homogeneous  over  the  array 
(due  to  the  array  designers  placing  the  plume  sources  near  the  beginnings  of  the  arrays, 
instead  of  the  middle)  or  that  the  model  is  missing  some  key  behaviour.  A  combination 
of  these  causes  is  possible. 

However,  that  the  model  can  describe  the  plume  well  without  an  array,  or  with  an 
array  with  individual  fits  for  different  downstream  distances  shows  it  captures  much  of  the 
key  physics. 

The  framework  developed  allows  a  link  between  previously  published  DSTO  work  on 

19This  problem  is  compounded  by  the  rescaling  of  source  concentrations  and  output  data  done  for 
long  upstream  distance  releases  (done  to  improve  detectability)  combined  with  an  unknown  algorithm  for 
cutting  off  noise  which  appears  only  partly  dependent  on  the  local  average  concentration. 
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fluctuating  plume  models  in  urban  settings  [30,  33],  current  DSTO-CSIRO  collaborative 
work  on  internal  concentration  fluctuations  [37],  and  more  complex  distributed  drag  CFD 
models  in  urban  areas  [18,  19,  20].  This  will  enable  a  high  fidelity  model  of  urban  dispersion 
to  be  implemented  for  emergency  response  applications,  with  the  extra  capability  for 
concentration  fluctuation  prediction  and  characterisation  of  the  uncertainty  in  modelling 
predictions.  Components  of  this  technology  may  ultimately  be  transferred  and  integrated 
into  a  national  Australian  hazard  prediction  modelling  system,  should  support  be  provided 
by  civilian  Government  agencies.20 

The  framework  developed  has  also  provided  a  link  through  the  UDM  to  another  impor¬ 
tant  operations  analysis  capability  that  DSTO  now  possesses — the  CBR  Virtual  Battle- 
space  developed  by  Dstl.  There  is  further  scope  for  development  of  concentration  reali¬ 
sation  models  within  this  capability  framework.  Such  models  will  rely  more  on  simpler, 
faster  approaches  to  dispersion  modelling  than  CFD  to  enable  rapid  calculations  for  oper¬ 
ational  analysis  tasks  relying  on  exploration  of  high  dimensional  variable  spaces,  thus  the 
framework  developed  here  is  ideal  for  these  purposes. 

Finally,  the  generation  of  prototype  simple  synthetic  environments  for  plume  realisa¬ 
tions  is  also  required  to  further  develop  data  fusion  algorithms  for  CBR  source  charac¬ 
terisation  and  network  detection.  The  framework  developed  here  is  the  first  step  in  the 
development  of  such  a  synthetic  environment  at  DSTO  for  these  purposes,  and  can  ulti¬ 
mately  be  used  to  compare  the  approaches  studied  at  DSTO  with  those  of  our  international 
collaborators. 


20A  project,  funded  by  Emergency  Management  Australia,  to  undertake  a  scoping  study  for  such  a 
capability  is  currently  underway  . 
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Appendix  A  New  Plume  Dispersion  Model 

A.l  Velocity  Profile 

The  parameter  a,  while  not  fundamental  (it  is  totally  determined  by  e  and  7)  is  important, 
because  the  formulae  for  U,  K,  d  and  20  cannot  be  written  in  terms  of  e.  Thus  we  have 
included  a  plot  of  its  behaviour  here  (see  Figure  Al).  The  key  features  it  has  are  that  it 
decreases  from  a  maximum  in  the  sparse  limit  and  approaches  0  at  infinity.  It  is  also  very 
sensitive  to  7,  increasing  sharply  as  7  increases. 

Following  on  from  this,  it  is  worthwhile  looking  in  more  detail  at  what  happens  to  the 
wind  profile  model  in  the  extremes  of  very  dense  and  very  sparse  canopies.  For  very  dense 
canopies  we  get  for  e  — >  00 
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7  ku*(z  —  H ). 

That  U+(z)  — >  00  in  this  limit  shows  that  this  theory  does  not  well  cover  the  case  of  an 
infinitely  dense  canopy.  This  is  an  artifact  of  the  functional  forms  forcing  the  velocity 
to  remain  constant  and  non-zero  at  H,  even  as  the  effective  ground  height  d  rises  to  H. 
To  negate  this,  it  is  necessary  to  let  the  velocity  at  the  canopy  top  decrease  to  zero  as  e 
increases  to  infinity,  such  that  the  velocity  at  the  top  of  the  ABL  remains  constant. 

In  the  limit  t->0a  different  problem  arises 
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The  problem  is  that  if  e  =  0  there  should  be  no  canopy,  thus  no  behaviour  change  in 
U{z )  at  z  =  H.  For  example  this  theory  predicts  a  significant  difference  between  wind 
passing  over  a  perfectly  flat  field  and  a  perfectly  flat  field  with  a  few  infinitesimally  thin 
sticks  standing  upon  it.  The  theory  behaves  sensibly  in  the  e  — >  0  limit  only  if  the  canopy 
height  H  — ►  0  at  the  same  time.  However  away  from  these  extremes,  the  U(z)  and  K(z) 
profiles  are  very  sensible  and  agree  well  with  experiment  (see  Section  4.4).  The  lower  limit 
of  e  for  which  we  can  use  this  theory  without  concern  is  determined  by  when  d  =  0.  The 
corresponding  e  value  is 

€  =  1.8336(7«)2.  (A3) 


A. 2  Vertical  concentration  gradient  within  the  canopy 


This  section  provides  a  closer  look  at  the  predicted  concentration  within  the  canopy. 
Inserting  the  modelled  formulae  for  U(z)  (41)  and  Kzz(z )  (38)  into  the  advection-diffusion 
equation  (21)  yields 


dC  _  <70  dC  d2C 

dx  2tanh {(3z/H)  dz  dz 2  ’ 


(A4) 


This  is  subject  to  the  boundary  condition  that  tracer  material  should  not  move  below  the 
ground: 

dC 


dz 


z= o 


For  typical  values  of  the  canopy  parameter  e  used  in  the  Coanda  experiment  (designed 
to  model  an  urban  environment),  tanh(/3z/Fl)  can  be  approximated  by  fdz/H.  This  is 
especially  true  near  the  ground.  The  solution  for  (A4)  can  thus  be  approximated  by 


C-(x,z)  ocx"3/4 
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exp 
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0  <  z  <  H. 


(A5) 


The  important  observations  to  make  here  are  that: 

A)  For  x  >  H/a  the  exponent  of  e  above  would  be  small  and  thus  C  within  the  canopy 
would  be  close  to  constant  a  short  distance  downstream  of  the  source, 

B)  The  exponent  a  (50)  should  indeed  be  2  at  the  ground. 
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cr 


Figure  Al:  The  parameter  a  against  e  for  various  7 
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